Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TRAIL.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 <fstream>
25 #include <string>
26 #include <list>
27 #include <map>
28 #include <sstream>
29 #include <iomanip>
30 #include <vector>
31 #include <cstdio>
32 #include <cmath>
33 #include <cstring>
34 #include <algorithm>
35 
36 #include "TRAIL_UnixUtils.H"
37 
38 //#ifdef _TRAIL_MPI_
39 #include "mpi.h"
40 //#else
41 //typedef int MPI_Comm;
42 //#endif
43 
44 #include "GEM.H"
45 #include "TRAIL.H"
46 
47 //#ifdef _ROCSTAR_X_
48 #include "roccom.h"
49 //#endif
50 
51 COM_EXTERN_MODULE( Rocface);
52 COM_EXTERN_MODULE( Rocout);
53 COM_EXTERN_MODULE( Rocin);
54 COM_EXTERN_MODULE( Rocblas);
55 COM_EXTERN_MODULE( Rocsurf);
56 COM_EXTERN_MODULE( Rocmop);
57 COM_EXTERN_MODULE( Rocprop);
58 
59 void
60 TRAIL_Debug(GEM_Partition &gp)
61 {
62  int rank = 0;
63  //#ifdef _TRAIL_MPI_
64  MPI_Comm_rank(gp._comm,&rank);
65  //#endif
66  if(!rank)
67  TRAIL_CreateDirectory("Roctrail");
68  //#ifdef _TRAIL_MPI_
69  MPI_Barrier(gp._comm);
70  //#endif
71  gp.debug(true);
72  std::ostringstream Ostr;
73  Ostr << "Roctrail/Roctrail_debug_" << gp._id;
74  std::ofstream *Ouf;
75  Ouf = new std::ofstream;
76  Ouf->open(Ostr.str().c_str());
77  gp._out = Ouf;
78 }
79 
80 // Get a string that encodes the given time. If the string is
81 // xx.yyyyyy, it means 0.yyyyyy*10^xx nanoseconds. This format
82 // ensures that sorting the strings alphabetically gives the right
83 // ordering of time. The string is null-terminated.
84 std::string
85 TRAIL_TimeString( double t)
86 {
87  std::ostringstream Ostr;
88  Ostr << std::scientific << std::setprecision(5) << t*1e10;
89  std::string buf(Ostr.str());
90  std::string timestring(buf.substr(9)+"."+buf.substr(0,1)+buf.substr(2,5));
91  return(timestring);
92 }
93 
94 double
95 TRAIL_TimeString(const std::string &ts)
96 {
97  std::string xp_s(ts.substr(0,2));
98  std::string num_s(ts.substr(3));
99  std::istringstream X_inst(ts);
100  int xp = 0;
101  X_inst >> xp;
102  xp = xp - 15;
103  double x = pow(10,(double)xp);
104  std::istringstream Num_inst(num_s);
105  double num = 0;
106  Num_inst >> num;
107  return(num*=x);
108 
109 }
110 
111 void
112 TRAIL_GetRocstarDumpStrings(const std::string &filename,std::string &wname,std::string &timestring,
113  std::string &rankstring)
114 {
115  std::string fname(filename.substr(0,filename.find_last_of(".")));
116  std::string::size_type x = fname.find_last_of("_");
117  std::string::size_type len = fname.length() - x;
118  rankstring.assign(fname.substr(x+1,len-1));
119  fname = fname.substr(0,x);
120  x = fname.find_last_of("_");
121  len = fname.length() - x;
122  timestring.assign(fname.substr(x+1,len-1));
123  wname = fname.substr(0,x);
124 }
125 
126 
127 //#ifdef _MESH_X_
128 int
129 TRAIL_ExtractPanes(const std::string &window_name,
130  const std::string &attribute_name,
131  int attribute_value,
132  std::vector<int> &pane_id)
133 {
134  std::vector<int> pane_ids;
135  pane_id.resize(0);
136  COM_get_panes(window_name.c_str(),pane_ids);
137  std::vector<int>::iterator pi = pane_ids.begin();
138  while(pi != pane_ids.end()){
139  int *attptr = NULL;
140  COM_get_array((window_name+"."+attribute_name).c_str(),*pi,&attptr);
141  if(*attptr == attribute_value)
142  pane_id.push_back(*pi);
143  pi++;
144  }
145  return 0;
146 }
147 
148 int
149 TRAIL_GetPanelAttribute(const std::string &window_name,
150  const std::string &attribute_name,
151  const std::string &qual_name,
152  int qualval,
153  std::vector<int> &attvec)
154 {
155  std::vector<int> pane_ids;
156  attvec.resize(0);
157  COM_get_panes(window_name.c_str(),pane_ids);
158  std::vector<int>::iterator pi = pane_ids.begin();
159  while(pi != pane_ids.end()){
160  int *attptr = NULL;
161  COM_get_array((window_name+"."+qual_name).c_str(),*pi,&attptr);
162  if(*attptr == qualval){
163  COM_get_array((window_name+"."+attribute_name).c_str(),*pi,&attptr);
164  if(attptr)
165  attvec.push_back(*attptr);
166  }
167  pi++;
168  }
169  return 0;
170 }
171 
172 int
173 TRAIL_GetPanelAttribute(const std::string &window_name,
174  const std::string &attribute_name,
175  std::vector<int> &attvec)
176 {
177  std::vector<int> pane_ids;
178  attvec.resize(0);
179  COM_get_panes(window_name.c_str(),pane_ids);
180  std::vector<int>::iterator pi = pane_ids.begin();
181  while(pi != pane_ids.end()){
182  int *attptr = NULL;
183  COM_get_array((window_name+"."+attribute_name).c_str(),*pi,&attptr);
184  if(attptr)
185  attvec.push_back(*attptr);
186  pi++;
187  }
188  return 0;
189 }
190 
191 int TRAIL_UniqueAcrossProcs(std::vector<int> &input_data,std::vector<int> &output_data,
192  MPI_Comm communicator)
193 {
194  int nprocs = 1;
195  MPI_Comm_size(communicator,&nprocs);
196  int nlocal_items = input_data.size();
197  std::vector<int> nitems_per_processor(nprocs);
198  MPI_Allgather(&nlocal_items,1,MPI_INT,&nitems_per_processor[0],1,MPI_INT,communicator);
199  std::vector<int> displacements(nprocs);
200  displacements[0] = 0;
201  int total_count = nitems_per_processor[0];
202  for(int i = 1;i < nprocs;i++){
203  total_count += nitems_per_processor[i];
204  displacements[i] = nitems_per_processor[i-1]+displacements[i-1];
205  }
206  std::vector<int> all_items(total_count);
207  MPI_Allgatherv(&input_data[0],nlocal_items,MPI_INT,
208  &all_items[0],&nitems_per_processor[0],&displacements[0],
209  MPI_INT,communicator);
210  std::sort(all_items.begin(),all_items.end());
211  std::vector<int>::iterator ui = std::unique(all_items.begin(),all_items.end());
212  output_data.resize(ui - all_items.begin());
213  std::copy(all_items.begin(),ui,output_data.begin());
214  return 0;
215 }
216 
228 int
229 TRAIL_Search_Block_Structured_Pool(std::vector<std::vector<int> > &search_extent,
230  std::vector<std::vector<std::vector<int> > > &extent_pool,
231  std::vector<std::vector<std::vector<int> > > &neighbor_extent,
232  std::vector<int> &neighbors)
233 {
234  size_t pool_depth = extent_pool.size();
235  size_t nd = search_extent.size();
236  neighbor_extent.resize(0);
237  neighbors.resize(0);
238  std::vector<std::vector<int> > extent(nd);
239  for(size_t i = 0;i < nd;i++)
240  extent[i].resize(2);
241  for(size_t i = 0;i < pool_depth;i++){
242  bool match = true;
243  for(size_t j = 0;j < nd;j++){
244  if(!((search_extent[j][0] >= extent_pool[i][j][0] &&
245  search_extent[j][0] <= extent_pool[i][j][1]) ||
246  (search_extent[j][1] >= extent_pool[i][j][0] &&
247  search_extent[j][1] <= extent_pool[i][j][1])))
248  match = false;
249  }
250  if(match){ // then partition has some searched for mesh points
251  neighbors.push_back(i);
252  for(size_t j = 0;j < nd;j++){
253  extent[j][0] = std::max(search_extent[j][0],extent_pool[i][j][0]);
254  extent[j][1] = std::min(search_extent[j][1],extent_pool[i][j][1]);
255  }
256  neighbor_extent.push_back(extent);
257  }
258  }
259  return(0);
260 }
261 
262 
276 int
277 TRAIL_Get_Block_Structured_Neighbors(std::vector<std::vector<int> > &local_extent,
278  std::vector<std::vector<int> > &global_extent,
279  std::vector<std::vector<std::vector<int> > > &extent_pool,
280  std::vector<std::vector<int> > &ghost_extent,
281  std::vector<std::vector<std::vector<int> > > &neighbor_extent,
282  std::vector<int> &neighbors)
283 {
284  size_t pool_depth = extent_pool.size();
285  size_t nd = local_extent.size();
286  std::vector<std::vector<std::vector<int> > > pool_extents(pool_depth);
287  neighbors.resize(0);
288  neighbor_extent.resize(0);
289  for(size_t i = 0;i < pool_depth;i++){
290  pool_extents[i].resize(nd);
291  for(size_t j = 0;j < nd;j++)
292  pool_extents[i][j].resize(2,0);
293  }
294  ghost_extent = local_extent;
295  for(size_t i = 0;i < nd;i++){
296  // Left neighbors
297  if(local_extent[i][0] > global_extent[i][0])
298  ghost_extent[i][0] = local_extent[i][0] - 1;
299  if(local_extent[i][1] < global_extent[i][1])
300  ghost_extent[i][1] = local_extent[i][1] + 1;
301  }
302  for(size_t i = 0;i < nd;i++){
303  // Left neighbors
304  if(local_extent[i][0] > global_extent[i][0]){
305  std::vector<std::vector<int> > search_extent(ghost_extent);
306  search_extent[i][1] = search_extent[i][0];
307  // search_extent[i][0] = search_extent[i][1] = ghost_extent[i][0];
308  std::vector<std::vector<std::vector<int> > > directional_extent;
309  std::vector<int> directional_neighbors;
310  TRAIL_Search_Block_Structured_Pool(search_extent,extent_pool,
311  directional_extent,directional_neighbors);
312  std::vector<int>::iterator dni = directional_neighbors.begin();
313  size_t ncount = 0;
314  while(dni != directional_neighbors.end()){
315  int neighbor_index = *dni++;
316  for(size_t j = 0;j < nd;j++){
317  if(pool_extents[neighbor_index][j][0] > 0) // This is done to avoid clash with previous setting
318  pool_extents[neighbor_index][j][0] = std::min(directional_extent[ncount][j][0],
319  pool_extents[neighbor_index][j][0]);
320  else
321  pool_extents[neighbor_index][j][0] = directional_extent[ncount][j][0];
322  if(pool_extents[neighbor_index][j][1] > 0)
323  pool_extents[neighbor_index][j][1] = std::max(directional_extent[ncount][j][1],
324  pool_extents[neighbor_index][j][1]);
325  else
326  pool_extents[neighbor_index][j][1] = directional_extent[ncount][j][1];
327 
328  }
329  ncount++;
330  // neighbors.push_back(neighbor_index);
331  }
332  }
333  // Right neighbors
334  if(local_extent[i][1] < global_extent[i][1]){
335  std::vector<std::vector<int> > search_extent(ghost_extent);
336  search_extent[i][0] = search_extent[i][1] = ghost_extent[i][1];
337  std::vector<std::vector<std::vector<int> > > directional_extent;
338  std::vector<int> directional_neighbors;
339  TRAIL_Search_Block_Structured_Pool(search_extent,extent_pool,
340  directional_extent,directional_neighbors);
341  std::vector<int>::iterator dni = directional_neighbors.begin();
342  size_t ncount = 0;
343  while(dni != directional_neighbors.end()){
344  int neighbor_index = *dni++;
345  for(size_t j = 0;j < nd;j++){
346  if(pool_extents[neighbor_index][j][0] > 0)
347  pool_extents[neighbor_index][j][0] = std::min(directional_extent[ncount][j][0],
348  pool_extents[neighbor_index][j][0]);
349  else
350  pool_extents[neighbor_index][j][0] = directional_extent[ncount][j][0];
351  if(pool_extents[neighbor_index][j][1] > 0)
352  pool_extents[neighbor_index][j][1] = std::max(directional_extent[ncount][j][1],
353  pool_extents[neighbor_index][j][1]);
354  else
355  pool_extents[neighbor_index][j][1] = directional_extent[ncount][j][1];
356 
357  }
358  ncount++;
359  // neighbors.push_back(neighbor_index);
360  }
361  }
362  }
363  for(size_t i = 0;i < pool_depth;i++){
364  if(pool_extents[i][0][0] > 0){
365  neighbors.push_back(i);
366  neighbor_extent.push_back(pool_extents[i]);
367  }
368  }
369  return(0);
370 }
371 
372 template<typename T>
373 void TRAIL_Copy2Attribute(const std::string &aname,const std::vector<T> &container,int pane_id,int asize = 1)
374 {
375 
376  COM_set_size(aname.c_str(),pane_id,asize);
377  COM_resize_array(aname.c_str(),pane_id);
378  T *leptr;
379  COM_get_array(aname.c_str(),pane_id,&leptr);
380  typename std::vector<T>::const_iterator ci = container.begin();
381  while(ci != container.end())
382  *leptr++ = *ci++;
383 }
384 
385 template<typename T>
386 void TRAIL_SetAttribute(const std::string &aname,int pane_id,T &value)
387 {
388 
389  COM_set_size(aname.c_str(),pane_id,1);
390  COM_resize_array(aname.c_str(),pane_id);
391  T *leptr;
392  COM_get_array(aname.c_str(),pane_id,&leptr);
393  *leptr = value;
394 }
395 
396 template<typename T>
397 void TRAIL_Copy2Attribute(const std::string &aname,const std::vector<std::vector<T> > &container,int pane_id)
398 {
399  unsigned int asize = container.size();
400  COM_set_size(aname.c_str(),pane_id,asize);
401  COM_resize_array(aname.c_str(),pane_id);
402  T *leptr;
403  COM_get_array(aname.c_str(),pane_id,&leptr);
404  typename std::vector<std::vector<T> >::const_iterator ci = container.begin();
405  while(ci != container.end()){
406  typename std::vector<T>::const_iterator vi = ci->begin();
407  while(vi != ci->end())
408  *leptr++ = *vi++;
409  ci++;
410  }
411 }
412 
414 int TRAIL_SurfaceMesh2Window(const std::string &wname,int pane_id,Mesh::NodalCoordinates &coords,
415  Mesh::Connectivity &conn)
416 {
417  COM_new_window(wname);
418  unsigned int number_of_nodes = coords.Size();
419  COM_set_size((wname+".nc").c_str(),pane_id,number_of_nodes);
420  COM_resize_array((wname+".nc").c_str(),pane_id);
421  double *nc = NULL;
422  COM_get_array((wname+".nc").c_str(),pane_id,&nc);
423  std::memcpy(nc,coords[1],number_of_nodes*3*sizeof(double));
424  Mesh::Connectivity tris;
425  Mesh::Connectivity quads;
426  Mesh::Connectivity::iterator ci = conn.begin();
427  bool known_element_type = false;
428  while(ci != conn.end()){
429  if(ci->size() == 3)
430  tris.AddElement(*ci++);
431  else if (ci->size() == 4)
432  quads.AddElement(*ci++);
433  else
434  assert(known_element_type);
435  }
436  unsigned int number_of_tris = tris.Nelem();
437  unsigned int number_of_quads = quads.Nelem();
438  int *con_array = NULL;
439  if(number_of_tris > 0){
440  COM_set_size((wname+".:t3:real").c_str(),pane_id,number_of_tris);
441  COM_resize_array((wname+".:t3:real").c_str(),pane_id);
442  COM_get_array((wname+".:t3:real").c_str(),pane_id,&con_array);
443  ci = tris.begin();
444  while(ci != tris.end()){
445  unsigned int element_id = ci - tris.begin() + 1;
446  unsigned int element_index = element_id - 1;
447  con_array[3*element_index] = (*ci)[0];
448  con_array[3*element_index+1] = (*ci)[1];
449  con_array[3*element_index+2] = (*ci)[2];
450  ci++;
451  }
452  }
453  if(number_of_quads > 0){
454  COM_set_size((wname+".:q4:").c_str(),pane_id,number_of_quads);
455  COM_resize_array((wname+".:q4:").c_str(),pane_id);
456  COM_get_array((wname+".:q4:").c_str(),pane_id,&con_array);
457  ci = quads.begin();
458  while(ci != quads.end()){
459  unsigned int element_id = ci - quads.begin() + 1;
460  unsigned int element_index = element_id - 1;
461  con_array[4*element_index] = (*ci)[0];
462  con_array[4*element_index+1] = (*ci)[1];
463  con_array[4*element_index+2] = (*ci)[2];
464  con_array[4*element_index+3] = (*ci)[3];
465  ci++;
466  }
467  }
468  COM_window_init_done(wname.c_str());
469  return(0);
470 }
471 
473 int TRAIL_UnstructuredMesh2Pane(const std::string &wname,int pane_id,
474  Mesh::UnstructuredMesh &mesh,
475  SolnMetaData &smdv,
476  std::vector<std::vector<double> > &soln_data,
477  int verblevel)
478 {
479 
480 
481 /*
482  * Note: This throws an error that the window is duplicated. Looking
483  * at Roccom_base::get_status it looks like
484  * Now, the window should be created prior to calling this function.
485  * George Zagaris (gzagaris@illinois.edu)
486  */
487 // if(COM_get_status(wname.c_str(),0)<0)
488 // COM_new_window(wname);
489 
490 
491  unsigned int number_of_nodes = mesh.nc.Size();
492 // std::cout << "Number of nodes: " << mesh.nc.Size() << std::endl;
493  COM_set_size((wname+".nc").c_str(),pane_id,number_of_nodes);
494  COM_resize_array((wname+".nc").c_str(),pane_id);
495  double *nc = NULL;
496  COM_get_array((wname+".nc").c_str(),pane_id,&nc);
497  std::memcpy(nc,mesh.nc[1],number_of_nodes*3*sizeof(double));
498  Mesh::Connectivity tets;
499  Mesh::Connectivity bricks;
500  Mesh::Connectivity::iterator ci = mesh.con.begin();
501 // std::cout << "Number of elements: " << mesh.con.Nelem( ) << std::endl;
502  bool known_element_type = false;
503  while(ci != mesh.con.end()){
504  if(ci->size() == 4)
505  tets.AddElement(*ci++);
506  else if (ci->size() == 8)
507  bricks.AddElement(*ci++);
508  else
509  assert(known_element_type);
510  }
511  unsigned int number_of_tets = tets.Nelem();
512  unsigned int number_of_bricks = bricks.Nelem();
513 // std::cout << "Number of tets: " << number_of_tets << std::endl;
514 // std::cout << "Number of bricks: " << number_of_bricks << std::endl;
515  int *con_array = NULL;
516  if(number_of_tets > 0){
517  COM_set_size((wname+".:T4:real").c_str(),pane_id,number_of_tets);
518  COM_resize_array((wname+".:T4:real").c_str(),pane_id);
519  COM_get_array((wname+".:T4:real").c_str(),pane_id,&con_array);
520  ci = tets.begin();
521  while(ci != tets.end()){
522  unsigned int element_id = ci - tets.begin() + 1;
523  unsigned int element_index = element_id - 1;
524  con_array[4*element_index] = (*ci)[0];
525  con_array[4*element_index+1] = (*ci)[1];
526  con_array[4*element_index+2] = (*ci)[2];
527  con_array[4*element_index+3] = (*ci)[3];
528  ci++;
529  }
530  }
531  if(number_of_bricks > 0){
532  COM_set_size((wname+".:B8:real").c_str(),pane_id,number_of_bricks);
533  COM_resize_array((wname+".:B8:real").c_str(),pane_id);
534  COM_get_array((wname+".:B8:real").c_str(),pane_id,&con_array);
535  ci = bricks.begin();
536  while(ci != bricks.end()){
537  unsigned int element_id = ci - bricks.begin() + 1;
538  unsigned int element_index = element_id - 1;
539  con_array[8*element_index] = (*ci)[0];
540  con_array[8*element_index+1] = (*ci)[1];
541  con_array[8*element_index+2] = (*ci)[2];
542  con_array[8*element_index+3] = (*ci)[3];
543  con_array[8*element_index+4] = (*ci)[4];
544  con_array[8*element_index+5] = (*ci)[5];
545  con_array[8*element_index+6] = (*ci)[6];
546  con_array[8*element_index+7] = (*ci)[7];
547  ci++;
548  }
549  }
550 
551  // Window is created and initialized outside
552  // COM_window_init_done(wname.c_str());
553 
554  return(0);
555 }
556 
564 int
565 TRAIL_FD2FE_WinCreate(const std::string &wname,const std::string &outwname,std::ostream *ouf)
566 {
567  // Get the window communicator
568  MPI_Comm communicator = MPI_COMM_NULL;
569  COM_get_communicator(wname.c_str(),&communicator);
570  IRAD::Comm::CommunicatorObject BaseComm(communicator);
571  int nprocs = BaseComm.Size();
572  int rank = BaseComm.Rank();
573  if(ouf && !rank)
574  *ouf << "TRAIL_AddBlockStructuredGhostZone::Nprocs = " << nprocs << std::endl;
575  std::vector<int> pane_ids;
576  COM_get_panes(wname.c_str(),pane_ids);
577  if(ouf && !rank)
578  *ouf << "Number of panes: " << pane_ids.size() << std::endl;
579  // Form a list of unique block id's across all processors
580  std::vector<int> local_block_ids;
581  TRAIL_GetPanelAttribute(wname,"block_id",local_block_ids);
582  if(ouf && !rank){
583  *ouf << "Local Block_ids: ";
584  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
585  *ouf << std::endl;
586  }
587  std::vector<int> global_blocks;
588  TRAIL_UniqueAcrossProcs(local_block_ids,global_blocks,BaseComm.World());
589  if(ouf && !rank){
590  *ouf << "Global Block_ids: ";
591  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
592  *ouf << std::endl;
593  }
594  // Now all_block_ids is an array of all the unique block ids across all procs
595  // For each block
596  std::vector<int>::iterator bii = global_blocks.begin();
597  std::string wname2 = outwname; // wname+"_uns";
598  bool window_created = false;
599  if(!window_created){
600  COM_new_window(wname2.c_str());
601  COM_new_attribute((wname2+".local_extent").c_str(), 'p',COM_INTEGER,6,"");
602  COM_new_attribute((wname2+".global_extent").c_str(),'p',COM_INTEGER,6,"");
603  COM_new_attribute((wname2+".send_extent").c_str(), 'p',COM_INTEGER,6,"");
604  COM_new_attribute((wname2+".send_panes").c_str(), 'p',COM_INTEGER,1,"");
605  COM_new_attribute((wname2+".recv_panes").c_str(), 'p',COM_INTEGER,1,"");
606  COM_new_attribute((wname2+".recv_extent").c_str(), 'p',COM_INTEGER,6,"");
607  // COM_new_attribute((wname2+".sister_pane").c_str(), 'p',COM_INTEGER,1,"");
608  COM_new_attribute((wname2+".block_id").c_str(), 'p',COM_INTEGER,1,"");
609  COM_new_attribute((wname2+".patch_id").c_str(), 'p',COM_INTEGER,1,"");
610  window_created = true;
611  COM_new_attribute((wname+".send_extent").c_str(), 'p',COM_INTEGER,6,"");
612  COM_new_attribute((wname+".send_panes").c_str(), 'p',COM_INTEGER,1,"");
613  COM_new_attribute((wname+".recv_extent").c_str(), 'p',COM_INTEGER,6,"");
614  COM_new_attribute((wname+".recv_panes").c_str(), 'p',COM_INTEGER,1,"");
615  // COM_new_attribute((wname+".sister_pane").c_str(), 'p',COM_INTEGER,1,"");
616  }
617  while(bii != global_blocks.end()){
618  int block_id = *bii++;
619  std::vector<int>::iterator fi = std::find(local_block_ids.begin(),local_block_ids.end(),block_id);
620  int block_color = 0;
621  if(fi != local_block_ids.end()) // Then this processor has data for the given block
622  block_color = 1;
623  // Split the communicator into haves and have nots for this block
624  if(ouf && !rank)
625  *ouf << "Splitting communicator for blocks." << std::endl;
626  IRAD::Comm::CommunicatorObject BlockComm;
627  BaseComm.Split(block_color,rank,BlockComm);
628  if(block_color == 1){ // all the guys with data on this block
629  int block_nproc = BlockComm.Size();
630  int block_rank = BlockComm.Rank();
631  // Form a list of unique patch id's across all processors
632  std::vector<int> local_patch_ids;
633  std::vector<int> panes;
634  TRAIL_ExtractPanes(wname,"block_id",block_id,panes); // get pane id's of panes with block id = block_id
635  std::vector<int>::iterator pi = panes.begin();
636  while(pi != panes.end()){
637  int *patch_id;
638  COM_get_array((wname+".patch_id").c_str(),*pi,&patch_id);
639  if(patch_id)
640  local_patch_ids.push_back(*patch_id);
641  pi++;
642  }
643  std::vector<int> all_patches(block_nproc);
644  TRAIL_UniqueAcrossProcs(local_patch_ids,all_patches,BlockComm.World());
645  // Now all_patches is an array of all the unique patch id across all procs having the current block
646  std::vector<int>::iterator pii = all_patches.begin();
647  while(pii != all_patches.end()){ // For each patch
648  int patch_id = *pii++;
649  // Determine if local processor owns part of the patch
650  // Split the communicator into haves and have nots
651  int patch_color = 0;
652  IRAD::Comm::CommunicatorObject PatchComm;
653  std::vector<int>::iterator fp = std::find(local_patch_ids.begin(),local_patch_ids.end(),patch_id);
654  if(fp != local_patch_ids.end())
655  patch_color = 1;
656  if(ouf && !rank)
657  *ouf << "Splitting communicator for patches." << std::endl;
658  BlockComm.Split(patch_color,block_rank,PatchComm);
659  if(patch_color == 1) { // all of us that have data on the given block/patch
660  int patch_nproc = PatchComm.Size();
661  int patch_rank = PatchComm.Rank();
662  std::vector<int> patch_pane;
663  TRAIL_ExtractPanes(wname,"patch_id",patch_id,patch_pane); // get the pane for this patch (hopefully only 1)
664  int *global_patch_extent = NULL;
665  COM_get_array((wname+".global_extent").c_str(),patch_pane[0],&global_patch_extent);
666  if(!global_patch_extent){
667  std::cerr << "ERROR: Window " << wname << " has no global_extent attribute." << std::endl;
668  return 1;
669  }
670  int *local_patch_extent_ptr = NULL;
671  COM_get_array((wname+".local_extent").c_str(),patch_pane[0],&local_patch_extent_ptr);
672  if(!local_patch_extent_ptr){
673  std::cerr << "ERROR: Window " << wname << " has no local_extent attribute." << std::endl;
674  return 1;
675  }
676  std::vector<int> local_patch_extent;
677  for(unsigned int aa = 0;aa < 6;aa++)
678  local_patch_extent.push_back(local_patch_extent_ptr[aa]);
679  if(ouf && !rank)
680  *ouf << "Getting Local Coordinates." << std::endl;
681  double *LocalCoordinates = NULL;
682  COM_get_array((wname+".nc").c_str(),patch_pane[0],&LocalCoordinates);
683  // Determine the total number of ghost nodes we need and augment the extended local patch extent
684  // and allocate an array for holding ALL the nodal coordinates
685  // Communicate the local patch extents
686  std::vector<int> all_patch_extents(6*patch_nproc);
687  PatchComm.AllGather(local_patch_extent,all_patch_extents,6,6);
688  std::vector<int> all_pane_ids(patch_nproc);
689  PatchComm.AllGather(patch_pane[0],all_pane_ids);
690  if(ouf && !rank)
691  *ouf << "All patch extents communicated." << std::endl;
692  std::vector<std::vector<std::vector<int> > > extent_pool(patch_nproc);
693  std::vector<std::vector<int> > ghost_extent;
694  std::vector<std::vector<std::vector<int> > > neighbor_extent;
695  std::vector<int> neighbors;
696  std::vector<std::vector<int> > local_extent(3);
697  std::vector<std::vector<int> > global_extent(3);
698  for(int dd = 0;dd < 3;dd++){
699  local_extent[dd].push_back(local_patch_extent[dd*2]);
700  local_extent[dd].push_back(local_patch_extent[dd*2+1]);
701  global_extent[dd].push_back(global_patch_extent[dd*2]);
702  global_extent[dd].push_back(global_patch_extent[dd*2+1]);
703  }
704  for(int pp = 0;pp < patch_nproc;pp++){
705  extent_pool[pp].resize(3);
706  for(int dd = 0;dd < 3;dd++){
707  extent_pool[pp][dd].push_back(all_patch_extents[pp*6+dd*2]);
708  extent_pool[pp][dd].push_back(all_patch_extents[pp*6+dd*2+1]);
709  }
710  }
711  if(ouf && !rank)
712  *ouf << "Calling TRAIL_Get_Block_Structured_Neighbors." << std::endl;
713  TRAIL_Get_Block_Structured_Neighbors(local_extent,
714  global_extent,
715  extent_pool,
716  ghost_extent,
717  neighbor_extent,
718  neighbors);
719 
720  if(ouf && !rank)
721  *ouf << "Calling TRAIL_Get_Block_Structured_Neighbors done." << std::endl;
722  // Loop thru neighbors and post receives for the Nodal Coordinates and
723  // the local extents requested from each neighbor.
724  std::vector<int>::iterator ni = neighbors.begin();
725  int number_of_neighbors = neighbors.size();
726  std::vector<std::vector<int> > RemoteNeighborExtent(number_of_neighbors);
727  std::vector<std::vector<double> > RemoteNodalCoordinates(number_of_neighbors);
728  std::vector<std::vector<double> > SendCoordinates(number_of_neighbors);
729  std::vector<std::vector<int> > FlattenedNeighborExtents(number_of_neighbors);
730  std::vector<int> RecvMsgID(number_of_neighbors);
731  PatchComm.Barrier();
732  if(ouf && !rank)
733  *ouf << "All processors Doing ComLoop 1" << std::endl;
734  for(int nc = 0;nc < number_of_neighbors;nc++){
735  int nid = *ni++;
736  if(nid != patch_rank){
737  RemoteNeighborExtent[nc].resize(6);
738  RecvMsgID[nc] = PatchComm.ARecv(RemoteNeighborExtent[nc],nid,1);
739  Mesh::BSExtent<int> bs_extent(neighbor_extent[nc]);
740  RemoteNodalCoordinates[nc].resize(3*bs_extent.NNodes());
741  bs_extent.Flatten(FlattenedNeighborExtents[nc]);
742  PatchComm.ASend(FlattenedNeighborExtents[nc],nid,1);
743  if(ouf){
744  *ouf << "Sending the extents I request from " << nid << ": ";
745  IRAD::Util::DumpContents(*ouf,FlattenedNeighborExtents[nc]," ");
746  *ouf << std::endl;
747  }
748  PatchComm.ARecv(RemoteNodalCoordinates[nc],nid,2);
749  }
750  }
751  // Loop thru the neighbors and send each one the indices we need and
752  // the nodal coordinates they need.
753  ni = neighbors.begin();
754  PatchComm.Barrier();
755  if(ouf && !rank)
756  *ouf << "Doing ComLoop 2." << std::endl;
757  for(int nc = 0;nc < number_of_neighbors;nc++){
758  int nid = *ni++;
759  if(nid != patch_rank){
760  PatchComm.WaitRecv(RecvMsgID[nc]);
761  Mesh::BSExtent<int> bs_extent(RemoteNeighborExtent[nc]);
762  int nnodes = bs_extent.NNodes();
763  if(ouf && !rank)
764  *ouf << "Sending " << nnodes << " nodal coordinates to remote proc." << std::endl;
765  SendCoordinates[nc].resize(0);
766  std::vector<int> flat_indices;
767  Mesh::BSExtent<int> LocalExtent(local_extent);
768  LocalExtent.GetFlatIndices(bs_extent,flat_indices);
769  std::vector<int>::iterator gii = flat_indices.begin();
770  while(gii != flat_indices.end()){
771  int nodeind = (*gii++ - 1)*3;
772  SendCoordinates[nc].push_back(LocalCoordinates[nodeind]);
773  SendCoordinates[nc].push_back(LocalCoordinates[nodeind + 1]);
774  SendCoordinates[nc].push_back(LocalCoordinates[nodeind + 2]);
775  }
776  if(ouf && !rank)
777  *ouf << "SendCoordinates.size == " << SendCoordinates[nc].size() << std::endl;
778  PatchComm.ASend(SendCoordinates[nc],nid,2);
779  }
780  }
781  PatchComm.Barrier();
782  if(ouf && !rank)
783  *ouf << "Waiting for all messages." << std::endl;
784  PatchComm.WaitAll();
785  PatchComm.ClearRequests();
786  if(ouf && !rank)
787  *ouf << "All messages received, requests cleared." << std::endl;
788  // Now all the coordinates we need to build the ghost zone are in RemoteNodalCoordinates array
789  // Here we register the new surface with Roccom. This will be the surface collecting points
790  // greedily on the right and not at all on the left.
791  // Create and size the unstructured extent array
792  Mesh::BSExtent<int> uns_extent(local_extent);
793  Mesh::BSExtent<int>::iterator uei = uns_extent.begin();
794  Mesh::BSExtent<int>::iterator gei = ghost_extent.begin();
795  while(uei != uns_extent.end())
796  (*uei++)[1] = (*gei++)[1];
797  uns_extent.Sync();
798 // if(ouf && !rank){
799 // std::vector<int> tempout;
800 // uns_extent.Flatten(tempout);
801 // *ouf << "Unstructured extent: ";
802 // IRAD::Util::DumpContents(*ouf,tempout," ");
803 // *ouf << std::endl;
804 // }
805  // Register it as the local_extent
806  std::vector<int> flatunsextent;
807  uns_extent.Flatten(flatunsextent);
808  TRAIL_Copy2Attribute((wname2+".local_extent"),flatunsextent,patch_pane[0]);
809  // Do the global
810  Mesh::BSExtent<int> GlobalExtent(global_extent);
811  std::vector<int> flatglobextent;
812  GlobalExtent.Flatten(flatglobextent);
813  TRAIL_Copy2Attribute((wname2+".global_extent"),flatunsextent,patch_pane[0]);
814  TRAIL_SetAttribute((wname2+".block_id"),patch_pane[0],block_id);
815  TRAIL_SetAttribute((wname2+".patch_id"),patch_pane[0],patch_id);
816  if(ouf && !rank)
817  *ouf << "Patch ID = " << patch_id << std::endl;
818  // TRAIL_SetAttribute((wname2+".sister_pane")),patch_pane[0],pane_id);
819  Mesh::Connectivity conn;
820  uns_extent.CreateUnstructuredMesh(conn);
821 // if(ouf && !rank)
822 // *ouf << "Root proc created " << conn.Nelem() << " elements." << std::endl
823 // << "Connectivity:" << std::endl
824 // << conn << std::endl;
825  std::vector<int> flatcon;
826  unsigned int nnodes_new = uns_extent.NNodes();
827  std::vector<double> NewCoords(3*nnodes_new);
828  conn.Flatten(flatcon);
829  COM_set_size((wname2+".nc").c_str(),patch_pane[0],nnodes_new);
830  COM_resize_array((wname2+".nc").c_str(),patch_pane[0]);
831  COM_set_size((wname2+".:q4:").c_str(),patch_pane[0],conn.Nelem());
832  COM_resize_array((wname2+".:q4:").c_str(),patch_pane[0]);
833  double *coords;
834  COM_get_array((wname2+".nc").c_str(),patch_pane[0],&coords);
835  int *conndat;
836  COM_get_array((wname2+".:q4:").c_str(),patch_pane[0],&conndat);
837  memcpy(conndat,&flatcon[0],sizeof(int)*4*conn.Nelem());
838  std::vector<int> localind;
839  Mesh::BSExtent<int> LocalExtent(local_extent);
840  unsigned int nlocalnodes = LocalExtent.NNodes();
841  uns_extent.GetFlatIndices(LocalExtent,localind);
842  for(unsigned int a = 0;a < nlocalnodes;a++){
843  coords[(localind[a]-1)*3] = LocalCoordinates[a*3];
844  coords[(localind[a]-1)*3+1] = LocalCoordinates[a*3+1];
845  coords[(localind[a]-1)*3+2] = LocalCoordinates[a*3+2];
846  }
847  std::vector<int> send_panes;
848  std::vector<int> recv_panes;
849  std::vector< std::vector<int> > flat_recv_extents;
850  std::vector< std::vector<int> > flat_send_extents;
851  for(int nc = 0;nc < number_of_neighbors;nc++){
852  int nid = neighbors[nc];
853  if(nid != patch_rank){
854  if(FlattenedNeighborExtents[nc][1] == LocalExtent[0][1]+1 ||
855  FlattenedNeighborExtents[nc][3] == LocalExtent[1][1]+1 ||
856  FlattenedNeighborExtents[nc][5] == LocalExtent[2][1]+1){
857  // Then it has coordinates on the right that I need
858  // First, trim the flattened neighbor extent so that it doesn't
859  // include any indices out of the local range
860  Mesh::BSExtent<int> NeighborExtent(FlattenedNeighborExtents[nc]);
861  Mesh::BSExtent<int> TrimmedNeighborExtent(NeighborExtent);
862  TrimmedNeighborExtent[0][0] = std::max(TrimmedNeighborExtent[0][0],uns_extent[0][0]);
863  TrimmedNeighborExtent[1][0] = std::max(TrimmedNeighborExtent[1][0],uns_extent[1][0]);
864  TrimmedNeighborExtent[2][0] = std::max(TrimmedNeighborExtent[2][0],uns_extent[2][0]);
865  std::vector<int> flatrecvextent;
866  TrimmedNeighborExtent.Flatten(flatrecvextent);
867  flat_recv_extents.push_back(flatrecvextent);
868  recv_panes.push_back(all_pane_ids[nid]);
869  // Now it's trimmed, and we want to copy the relavent coordinates from
870  // the original untrimmed received coordinates array to the appropriate
871  // spot in the new coordinates array. To do so, we need the flat indices
872  // of the TrimmedNeighborExtent wrt both the source and destination.
873  std::vector<int> source_indices;
874  std::vector<int> dest_indices;
875  NeighborExtent.GetFlatIndices(TrimmedNeighborExtent,source_indices);
876  uns_extent.GetFlatIndices(TrimmedNeighborExtent,dest_indices);
877  unsigned int ncopynodes = source_indices.size();
878  for(unsigned int a = 0;a < ncopynodes;a++){
879  unsigned int sind = (source_indices[a]-1)*3;
880  unsigned int dind = (dest_indices[a]-1)*3;
881  coords[dind] = RemoteNodalCoordinates[nc][sind];
882  coords[dind+1] = RemoteNodalCoordinates[nc][sind+1];
883  coords[dind+2] = RemoteNodalCoordinates[nc][sind+2];
884  }
885  }
886  if(FlattenedNeighborExtents[nc][1] == LocalExtent[0][0]-1 ||
887  FlattenedNeighborExtents[nc][3] == LocalExtent[1][0]-1 ||
888  FlattenedNeighborExtents[nc][5] == LocalExtent[2][0]-1){
889  if(ouf)
890  *ouf << "Yeah, I send left to pane " << all_pane_ids[nid] << std::endl;
891  // Then this processor is on the left
892  // First, trim the flattened neighbor extent so that it doesn't
893  // include any indices out of the local range
894  send_panes.push_back(all_pane_ids[nid]);
895  Mesh::BSExtent<int> NeighborExtent(RemoteNeighborExtent[nc]); // extent requested by remote proc
896  // Mesh::BSExtent<int> TrimmedNeighborExtent(NeighborExtent);
897  // TrimmedNeighborExtent[0][0] = std::min(TrimmedNeighborExtent[0][1],uns_extent[0][1]);
898  // TrimmedNeighborExtent[1][0] = std::min(TrimmedNeighborExtent[1][1],uns_extent[1][1]);
899  // TrimmedNeighborExtent[2][0] = std::min(TrimmedNeighborExtent[2][1],uns_extent[2][1]);
900  std::vector<int> flatsendextent;
901  NeighborExtent.Flatten(flatsendextent);
902  if(ouf){
903  *ouf << "Send extents: ";
904  IRAD::Util::DumpContents(*ouf,flatsendextent," ");
905  *ouf << std::endl;
906  }
907  flat_send_extents.push_back(flatsendextent);
908  }
909  }
910  }
911  int nrecv = recv_panes.size();
912  int nsend = send_panes.size();
913  TRAIL_Copy2Attribute((wname+".send_panes"),send_panes,patch_pane[0],nsend);
914  TRAIL_Copy2Attribute((wname2+".send_panes"),send_panes,patch_pane[0],nsend);
915  TRAIL_Copy2Attribute((wname+".recv_panes"),recv_panes,patch_pane[0],nrecv);
916  TRAIL_Copy2Attribute((wname2+".recv_panes"),recv_panes,patch_pane[0],nrecv);
917  TRAIL_Copy2Attribute((wname+".send_extent"),flat_send_extents,patch_pane[0]);
918  TRAIL_Copy2Attribute((wname2+".recv_extent"),flat_recv_extents,patch_pane[0]);
919  // COM_window_init_done(wname.c_str());
920  // COM_window_init_done(wname2.c_str());
921  if(ouf && !rank){
922  *ouf << "Number of panes I send to: " << nsend << std::endl
923  << "Number of panes I recv from: " << nrecv << std::endl
924  << "Send panes: ";
925  IRAD::Util::DumpContents(*ouf,send_panes," ");
926  *ouf << std::endl
927  << "Recv panes: ";
928  IRAD::Util::DumpContents(*ouf,recv_panes," ");
929  *ouf << std::endl
930  << "Send Extents: " << std::endl;
931  for(int o = 0;o < nsend;o++){
932  IRAD::Util::DumpContents(*ouf,flat_send_extents[o]," ");
933  *ouf << std::endl;
934  }
935  *ouf << "Recv Extents:" << std::endl;
936  for(int o = 0;o < nrecv;o++){
937  IRAD::Util::DumpContents(*ouf,flat_recv_extents[o]," ");
938  *ouf << std::endl;
939  }
940  }
941  }
942  }
943  // MPI_Free_Comm (if needed, happens automatically when CommunicatorObject is destroyed)
944  }
945  // MPI_Free_Comm (if needed, happens automatically when CommunicatorObject is destroyed)
946  }
947  BaseComm.Barrier();
948  return(0);
949 }
950 
951 template<typename T>
952 int TRAIL_Att2Vec(const std::string &att,int pane_id,std::vector<T> &dest)
953 {
954  T *ptr = NULL;
955  dest.resize(0);
956  int ncomp = 0;
957  std::string unit;
958  int insize = 0;
959  COM_Type type;
960  char loc;
961  COM_get_attribute(att.c_str(),&loc,&type,&ncomp,&unit);
962  COM_get_size(att.c_str(),pane_id,&insize);
963  COM_get_array(att.c_str(),pane_id,&ptr);
964  if(!ptr)
965  return 1;
966  int count = 0;
967  int total = ncomp*insize;
968  while(count < total)
969  dest.push_back(ptr[count++]);
970  return 0;
971 }
972 
973 int
974 TRAIL_FE2FD_Transfer(const std::string &fewin,const std::string &fdwin,
975  const std::string &attlist,MPI_Comm communicator,
976  std::ostream *ouf)
977 {
978  std::vector<std::string> atts;
979  std::istringstream Istr(attlist);
980  std::string attstr;
981  if(ouf)
982  *ouf << "FE2FD_Transfer: Enter" << std::endl;
983  MPI_Comm oldcomm = COM_get_default_communicator();
984  COM_set_default_communicator(communicator);
985  while(Istr >> attstr)
986  atts.push_back(attstr);
987  // Get the window communicator
988  // MPI_Comm communicator = MPI_COMM_NULL;
989  // COM_get_communicator(fdwin.c_str(),&communicator);
990  IRAD::Comm::CommunicatorObject BaseComm(communicator);
991  int rank = BaseComm.Rank();
992  std::vector<std::string>::iterator atti = atts.begin();
993  if(ouf)
994  *ouf << "FE2FD: Reducing attributes over shared nodes." << std::endl;
995  while(atti != atts.end()){
996  std::string attstr = fewin+"."+*atti++;
997  if(ouf)
998  *ouf << "FE2FD: Reducing " << attstr << std::endl;
999  int att_hndl = COM_get_attribute_handle(attstr);
1000  int att_color = 0;
1001  if(att_hndl >= 0)
1002  att_color = 1;
1003  IRAD::Comm::CommunicatorObject AttComm;
1004  BaseComm.Split(att_color,rank,AttComm);
1005  if(att_color == 1){
1006  MPI_Comm ocomm = COM_get_default_communicator();
1007  COM_set_default_communicator(AttComm.GetCommunicator());
1008  COM_LOAD_MODULE_STATIC_DYNAMIC(Rocmap,"MAP");
1009  int mean_over_shared_nodes = COM_get_function_handle("MAP.reduce_average_on_shared_nodes");
1010  COM_call_function(mean_over_shared_nodes,&att_hndl);
1011  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocmap,"MAP");
1012  COM_set_default_communicator(ocomm);
1013  if(ouf)
1014  *ouf << "FE2FD: Reduced " << attstr << std::endl;
1015  }
1016  else if (ouf)
1017  *ouf << "FE2FD: " << attstr << " did not exist." << std::endl;
1018  }
1019  // COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocmap,"MAP");
1020  if(ouf)
1021  *ouf << "FE2FD: Waiting on all procs to finish reducing." << std::endl;
1022  BaseComm.Barrier();
1023  if(ouf)
1024  *ouf << "FE2FD: All procs done reducing attributes." << std::endl;
1025  // int nprocs = BaseComm.Size();
1026  std::vector<int> pane_ids;
1027  COM_get_panes(fdwin.c_str(),pane_ids);
1028  if(ouf)
1029  *ouf << "FE2FD: Number of panes: " << pane_ids.size() << std::endl;
1030  // Form a list of unique block id's across all processors
1031  std::vector<int> local_block_ids;
1032  TRAIL_GetPanelAttribute(fdwin,"block_id",local_block_ids);
1033  if(ouf){
1034  *ouf << "FE2FD: Local Block_ids: ";
1035  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
1036  *ouf << std::endl;
1037  }
1038  std::vector<int> global_blocks;
1039  TRAIL_UniqueAcrossProcs(local_block_ids,global_blocks,BaseComm.World());
1040  if(ouf){
1041  *ouf << "FE2FD: Global Block_ids: ";
1042  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
1043  *ouf << std::endl;
1044  }
1045  // Now all_block_ids is an array of all the unique block ids across all procs
1046  // For each block
1047  std::vector<int>::iterator bii = global_blocks.begin();
1048  while(bii != global_blocks.end()){
1049  int block_id = *bii++;
1050  std::vector<int>::iterator fi = std::find(local_block_ids.begin(),local_block_ids.end(),block_id);
1051  int block_color = 0;
1052  if(fi != local_block_ids.end()) // Then this processor has data for the given block
1053  block_color = 1;
1054  // Split the communicator into haves and have nots for this block
1055  if(ouf)
1056  *ouf << "FE2FD: Splitting communicator for blocks." << std::endl;
1057  IRAD::Comm::CommunicatorObject BlockComm;
1058  BaseComm.Split(block_color,rank,BlockComm);
1059  if(block_color == 1){ // all the guys with data on this block
1060  int block_nproc = BlockComm.Size();
1061  int block_rank = BlockComm.Rank();
1062  // Form a list of unique patch id's across all processors
1063  std::vector<int> local_patch_ids;
1064  std::vector<int> panes;
1065  TRAIL_ExtractPanes(fdwin,"block_id",block_id,panes); // get pane id's of panes with block id = block_id
1066  std::vector<int>::iterator pi = panes.begin();
1067  while(pi != panes.end()){
1068  int *patch_id;
1069  COM_get_array((fdwin+".patch_id").c_str(),*pi,&patch_id);
1070  if(patch_id)
1071  local_patch_ids.push_back(*patch_id);
1072  pi++;
1073  }
1074  std::vector<int> all_patches(block_nproc);
1075  TRAIL_UniqueAcrossProcs(local_patch_ids,all_patches,BlockComm.World());
1076  // Now all_patches is an array of all the unique patch id across all procs having the current block
1077  std::vector<int>::iterator pii = all_patches.begin();
1078  while(pii != all_patches.end()){ // For each patch
1079  int patch_id = *pii++;
1080  // Determine if local processor owns part of the patch
1081  // Split the communicator into haves and have nots
1082  int patch_color = 0;
1083  IRAD::Comm::CommunicatorObject PatchComm;
1084  std::vector<int>::iterator fp = std::find(local_patch_ids.begin(),local_patch_ids.end(),patch_id);
1085  if(fp != local_patch_ids.end())
1086  patch_color = 1;
1087  if(ouf)
1088  *ouf << "FE2FD: Splitting communicator for patches." << std::endl;
1089  BlockComm.Split(patch_color,block_rank,PatchComm);
1090  if(patch_color == 1) { // all of us that have data on the given block/patch
1091  // int patch_nproc = PatchComm.Size();
1092  // int patch_rank = PatchComm.Rank();
1093  std::vector<int> patch_pane;
1094  TRAIL_ExtractPanes(fdwin,"patch_id",patch_id,patch_pane); // get the pane for this patch (hopefully only 1)
1095  // std::vector<int>::iterator pii = pane_ids.begin();
1096  // while(pii != pane_ids.end()){
1097  int pane_id = patch_pane[0];
1098  int err = 0;
1099  std::vector<int> flat_uns_extent;
1100  err = TRAIL_Att2Vec(fewin+".local_extent",pane_id,flat_uns_extent);
1101  std::vector<int> flat_fd_extent;
1102  err = TRAIL_Att2Vec(fdwin+".local_extent",pane_id,flat_fd_extent);
1103  assert(err == 0);
1104  Mesh::BSExtent<int> UnsExtent(flat_uns_extent);
1105  Mesh::BSExtent<int> FDExtent(flat_fd_extent);
1106  std::vector<int> src_indices;
1107  UnsExtent.GetFlatIndices(FDExtent,src_indices);
1108  std::vector<int> asizes(atts.size(),0);
1109  int attindex = 0;
1110  std::vector<std::string>::iterator atti = atts.begin();
1111  while(atti != atts.end()){
1112  std::string att = *atti++;
1113  int ncomp = 0;
1114  char loc;
1115  COM_Type type;
1116  std::string unit;
1117  COM_get_attribute((fewin+"."+att).c_str(),&loc,&type,&ncomp,&unit);
1118  asizes[attindex++] = COM_get_sizeof(type,ncomp);
1119  }
1120  if(ouf)
1121  *ouf << "FE2FD: All attributes sized up." << std::endl;
1122 
1123  atti = atts.begin();
1124  attindex = 0;
1125  if(ouf)
1126  *ouf << "FE2FD: Copying buffers." << std::endl;
1127  while(atti != atts.end()){
1128  std::string attstr = fewin+"."+*atti;
1129  std::string trgatt = fdwin+"."+*atti++;
1130  char *srcptr = NULL;
1131  char *trgptr = NULL;
1132  COM_get_array(attstr.c_str(),pane_id,&srcptr);
1133  COM_get_array(trgatt.c_str(),pane_id,&trgptr);
1134  std::vector<int>::iterator sii = src_indices.begin();
1135  int count = 0;
1136  while(sii != src_indices.end()){
1137  int src_index = *sii++ - 1;
1138  memcpy(&trgptr[count++*asizes[attindex]],&srcptr[src_index*asizes[attindex]],asizes[attindex]);
1139  }
1140  attindex++;
1141  }
1142  }
1143  }
1144  }
1145  }
1146  if(ouf)
1147  *ouf << "FD2FE: Waiting at final barrier" << std::endl;
1148  BaseComm.Barrier();
1149  if(ouf)
1150  *ouf << "FD2FE: Exit" << std::endl;
1151  COM_set_default_communicator(oldcomm);
1152  return(0);
1153 }
1154 
1155 int
1156 TRAIL_Add_Attributes(const std::string &srcwin,
1157  const std::string &destwin,
1158  const std::vector<std::string> &atts)
1159 {
1160  std::vector<std::string>::const_iterator ai = atts.begin();
1161  while(ai != atts.end()){
1162  std::string att = *ai++;
1163  std::string tatt = destwin+"."+att;
1164  std::string satt = srcwin+"."+att;
1165  // First, make sure the destination attribute exists, if not then create and size it
1166  {
1167  if(COM_get_attribute_handle(tatt) <= 0){
1168  int ncomp = 0;
1169  char loc;
1170  COM_Type type;
1171  std::string unit;
1172  COM_get_attribute(satt,&loc,&type,&ncomp,&unit);
1173  COM_new_attribute(tatt,loc,type,ncomp,unit.c_str());
1174  // COM_window_init_done(destwin,0);
1175  }
1176  }
1177  }
1178  COM_window_init_done(destwin,0);
1179  return(0);
1180 }
1181 
1182 int
1183 TRAIL_FD2FE_Transfer(const std::string &srcwin,const std::string &destwin,
1184  const std::string &attlist,std::ostream *ouf)
1185 {
1186  std::vector<std::string> atts;
1187  std::istringstream Istr(attlist);
1188  std::string attstr;
1189  while(Istr >> attstr)
1190  atts.push_back(attstr);
1191  // Get the window communicator
1192  MPI_Comm communicator = MPI_COMM_NULL;
1193  COM_get_communicator(srcwin.c_str(),&communicator);
1194  IRAD::Comm::CommunicatorObject BaseComm(communicator);
1195  IRAD::Comm::CommunicatorObject BlockComm;
1196  IRAD::Comm::CommunicatorObject PatchComm;
1197  // Make sure the attribute exists
1198  TRAIL_Add_Attributes(srcwin,destwin,atts);
1199  // int nprocs = BaseComm.Size();
1200  int rank = BaseComm.Rank();
1201  std::vector<int> pane_ids;
1202  COM_get_panes(srcwin.c_str(),pane_ids);
1203  if(ouf && !rank)
1204  *ouf << "TRAIL: Number of panes: " << pane_ids.size() << std::endl;
1205  // Form a list of unique block id's across all processors
1206  std::vector<int> local_block_ids;
1207  TRAIL_GetPanelAttribute(srcwin,"block_id",local_block_ids);
1208  if(ouf && !rank){
1209  *ouf << "TRAIL: Local Block_ids: ";
1210  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
1211  *ouf << std::endl;
1212  }
1213  std::vector<int> global_blocks;
1214  TRAIL_UniqueAcrossProcs(local_block_ids,global_blocks,BaseComm.World());
1215  if(ouf && !rank){
1216  *ouf << "TRAIL: Global Block_ids: ";
1217  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
1218  *ouf << std::endl;
1219  }
1220  // Now all_block_ids is an array of all the unique block ids across all procs
1221  // For each block
1222  if(true){
1223  std::vector<int>::iterator bii = global_blocks.begin();
1224  while(bii != global_blocks.end()){
1225  int block_id = *bii++;
1226  std::vector<int>::iterator fi = std::find(local_block_ids.begin(),local_block_ids.end(),block_id);
1227  int block_color = 0;
1228  if(fi != local_block_ids.end()) // Then this processor has data for the given block
1229  block_color = 1;
1230  // Split the communicator into haves and have nots for this block
1231  if(ouf && !rank)
1232  *ouf << "TRAIL: Splitting communicator for blocks." << std::endl;
1233  // IRAD::Comm::CommunicatorObject BlockComm;
1234  BaseComm.Split(block_color,rank,BlockComm);
1235  if(true){
1236  if(block_color == 1){ // all the guys with data on this block
1237  int block_nproc = BlockComm.Size();
1238  int block_rank = BlockComm.Rank();
1239  // Form a list of unique patch id's across all processors
1240  std::vector<int> local_patch_ids;
1241  std::vector<int> panes;
1242  TRAIL_ExtractPanes(srcwin,"block_id",block_id,panes); // get pane id's of panes with block id = block_id
1243  std::vector<int>::iterator pi = panes.begin();
1244  while(pi != panes.end()){
1245  int *patch_id;
1246  COM_get_array((srcwin+".patch_id").c_str(),*pi,&patch_id);
1247  if(patch_id)
1248  local_patch_ids.push_back(*patch_id);
1249  pi++;
1250  }
1251  std::vector<int> all_patches(block_nproc);
1252  TRAIL_UniqueAcrossProcs(local_patch_ids,all_patches,BlockComm.World());
1253  // Now all_patches is an array of all the unique patch id across all procs having the current block
1254  std::vector<int>::iterator pii = all_patches.begin();
1255  while(pii != all_patches.end()){ // For each patch
1256  int patch_id = *pii++;
1257  // Determine if local processor owns part of the patch
1258  // Split the communicator into haves and have nots
1259  int patch_color = 0;
1260  // IRAD::Comm::CommunicatorObject PatchComm;
1261  std::vector<int>::iterator fp = std::find(local_patch_ids.begin(),local_patch_ids.end(),patch_id);
1262  if(fp != local_patch_ids.end())
1263  patch_color = 1;
1264  if(ouf && !rank)
1265  *ouf << "TRAIL: Splitting communicator for patches." << std::endl;
1266  BlockComm.Split(patch_color,block_rank,PatchComm);
1267  if(true){
1268  if(patch_color == 1) { // all of us that have data on the given block/patch
1269  int patch_nproc = PatchComm.Size();
1270  // int patch_rank = PatchComm.Rank();
1271  std::vector<int> patch_pane;
1272  TRAIL_ExtractPanes(srcwin,"patch_id",patch_id,patch_pane); // get the pane for this patch (hopefully only 1)
1273  // std::vector<int>::iterator pii = pane_ids.begin();
1274  // while(pii != pane_ids.end()){
1275  int pane_id = patch_pane[0];
1276  std::vector<int> all_pane_ids(patch_nproc);
1277  PatchComm.AllGather(pane_id,all_pane_ids);
1278  int err = 0;
1279  std::vector<int> flat_uns_extent;
1280  err = TRAIL_Att2Vec(destwin+".local_extent",pane_id,flat_uns_extent);
1281  std::vector<int> flat_global_extent;
1282  err = TRAIL_Att2Vec(destwin+".global_extent",pane_id,flat_global_extent);
1283  std::vector<int> shared_panes;
1284  err = TRAIL_Att2Vec(destwin+".shared_panes",pane_id,shared_panes);
1285  std::vector<int> flat_shared_extents;
1286  err = TRAIL_Att2Vec(destwin+".shared_extents",pane_id,flat_shared_extents);
1287  std::vector<int> flat_fd_extent;
1288  err = TRAIL_Att2Vec(srcwin+".local_extent",pane_id,flat_fd_extent);
1289  assert(err == 0);
1290 
1291  if(ouf){
1292  *ouf << "fd_extent: (["
1293  << flat_fd_extent[0] << "," << flat_fd_extent[1] << "],["
1294  << flat_fd_extent[2] << "," << flat_fd_extent[3] << "],["
1295  << flat_fd_extent[4] << "," << flat_fd_extent[5] << "])"
1296  << std::endl;
1297  *ouf << "unstructured_extent: (["
1298  << flat_uns_extent[0] << "," << flat_uns_extent[1] << "],["
1299  << flat_uns_extent[2] << "," << flat_uns_extent[3] << "],["
1300  << flat_uns_extent[4] << "," << flat_uns_extent[5] << "])"
1301  << std::endl;
1302  *ouf << "global_extent: (["
1303  << flat_global_extent[0] << "," << flat_global_extent[1] << "],["
1304  << flat_global_extent[2] << "," << flat_global_extent[3] << "],["
1305  << flat_global_extent[4] << "," << flat_global_extent[5] << "])"
1306  << std::endl;
1307  for(int iii = 0;iii < flat_shared_extents.size()/6;iii++){
1308  unsigned int index = iii*6;
1309  *ouf << "shared_extents (" << shared_panes[iii] << "): (["
1310  << flat_shared_extents[index+0] << "," << flat_shared_extents[index+1] << "],["
1311  << flat_shared_extents[index+2] << "," << flat_shared_extents[index+3] << "],["
1312  << flat_shared_extents[index+4] << "," << flat_shared_extents[index+5] << "])"
1313  << std::endl;
1314  }
1315  }
1316 
1317  std::vector<int> all_fd_extents(6*patch_nproc);
1318  // PatchComm.AllGather<std::vector<int>,int>(flat_fd_extent,all_fd_extents,6,6);
1319  PatchComm.AllGather(flat_fd_extent,all_fd_extents,6,6);
1320  if(ouf && !rank)
1321  *ouf << "TRAIL: All FD Extents communicated." << std::endl;
1322  // Now need to determine the interpane communication lists so that
1323  // we can gather the FD representation data needed to populate the
1324  // full extent of the FE representation.
1325  Mesh::BSExtent<int> UnsExtent(flat_uns_extent);
1326  Mesh::BSExtent<int> GlobalExtent(flat_global_extent);
1327  Mesh::BSExtent<int> FDExtent(flat_fd_extent);
1328 
1329  // Get the indices of the original fd gridpoints wrt the fe mesh.
1330  std::vector<int> theflatindices;
1331  UnsExtent.GetFlatIndices(FDExtent,theflatindices);
1332 
1333  std::vector<Mesh::BSExtent<int> > AllFDExtents;
1334  std::vector<Mesh::BSExtent<int> > SharedExtents;
1335  std::vector<Mesh::BSExtent<int> > RecvExtents;
1336  std::vector<int> remote_ranks;
1337  // std::vector<int>::iterator spi = shared_panes.begin();
1338  std::vector<std::vector<char> > SendBuffers;
1339  std::vector<std::vector<char> > RecvBuffers;
1340  std::vector<std::vector<int> > FlatIndices;
1341  std::vector<int> recvranks;
1342  std::vector<int> sendranks;
1343  for(int i = 0;i < patch_nproc;i++)
1344  AllFDExtents.push_back(Mesh::BSExtent<int>(&all_fd_extents[i*6]));
1345 
1346  int count = 0;
1347  int nsend = 0;
1348  int nrecv = 0;
1349  std::vector<std::string>::iterator ai = atts.begin();
1350  std::vector<int> asizes(atts.size(),0);
1351  size_t bsize = 0; // block size = all attribute sizes added up
1352  int attindex = 0;
1353  while(ai != atts.end()){
1354  std::string att = *ai++;
1355  int ncomp = 0;
1356  char loc;
1357  COM_Type type;
1358  std::string unit;
1359  if(ouf && !rank)
1360  *ouf << "TRAIL: Processing attribute " << att << std::endl;
1361  COM_get_attribute((srcwin+"."+att).c_str(),&loc,&type,&ncomp,&unit);
1362  asizes[attindex] = COM_get_sizeof(type,ncomp);
1363  bsize += asizes[attindex++];
1364  }
1365  if(ouf && !rank)
1366  *ouf << "TRAIL: All attributes sized up." << std::endl;
1367  std::vector<int>::iterator spi = shared_panes.begin();
1368  while(spi != shared_panes.end()){
1369  Mesh::BSExtent<int> SharedExtent(&flat_shared_extents[6*count++]);
1370  SharedExtents.push_back(SharedExtent);
1371  int rpane_id = *spi++;
1372  // If the shared extent (FE rep) overlaps with
1373  // the FDExtent, then those nodes are send nodes
1374  int remote_rank = std::find(all_pane_ids.begin(),all_pane_ids.end(),rpane_id) - all_pane_ids.begin();
1375  remote_ranks.push_back(remote_rank);
1376  Mesh::BSExtent<int> CommExtent;
1377  FDExtent.Overlap(SharedExtent,CommExtent);
1378  if(!CommExtent.empty()){ // Then we need to send data to the remote processor
1379  // For each attribute to send, pack and send the buffer (which needs to be persistent)
1380  sendranks.push_back(remote_rank);
1381  int nnodes = CommExtent.NNodes();
1382  std::vector<std::string>::iterator ai = atts.begin();
1383  std::vector<int> indices;
1384  FDExtent.GetFlatIndices(CommExtent,indices);
1385  // SendBuffers.push_back(std::vector<char>(nnodes*bsize,0));
1386  std::vector<char> sendbufv(nnodes*bsize,0);
1387  SendBuffers.push_back(sendbufv);
1388  std::vector<int>::iterator ii = indices.begin();
1389  size_t byte_offset = 0;
1390  while(ii != indices.end()){
1391  int attindex = 0;
1392  int index = *ii++ - 1;
1393  ai = atts.begin();
1394  while(ai != atts.end()){
1395  std::string att = *ai++;
1396  char *data_ptr = NULL;
1397  COM_get_array((srcwin+"."+att).c_str(),pane_id,&data_ptr);
1398  // memcpy(&SendBuffers[nsend][byte_offset],&data_ptr[(index-1)*asizes[attindex]],asizes[attindex]);
1399  memcpy(&SendBuffers[nsend][byte_offset],&data_ptr[index*asizes[attindex]],asizes[attindex]);
1400 
1401  byte_offset += asizes[attindex++];
1402  }
1403  }
1404  // PatchComm.ASend(SendBuffers[nsend++],remote_rank);
1405  nsend++;
1406  }
1407  else { // Then we need to receive data from the remote processor
1408  // First, use the overlap function to ensure that we build a receive
1409  // buffer for only the *real* nodes on the remote processor. Each
1410  // processor sends only real nodes as seen in the above send block
1411  AllFDExtents[remote_rank].Overlap(SharedExtent,CommExtent);
1412  if(!CommExtent.empty()){ // if it's empty, then there were only ghosts, no receive necessary
1413  RecvExtents.push_back(CommExtent);
1414  // For each attribute to receive build a buffer into which to receive
1415  int nnodes = CommExtent.NNodes();
1416  std::vector<char> recvbuf(nnodes*bsize,0);
1417  std::vector<std::string>::iterator ai = atts.begin();
1418  // RecvBuffers.push_back(std::vector<char>(nnodes*bsize,0));
1419  RecvBuffers.push_back(recvbuf);
1420  if(ouf && !rank)
1421  *ouf << "TRAIL: Receive buffer size for " << nnodes << " nodes is "
1422  << nnodes*bsize << "." << std::endl;
1423  recvranks.push_back(remote_rank);
1424  // PatchComm.ARecv(RecvBuffers[nrecv++],remote_rank);
1425  }
1426  }
1427  }
1428  PatchComm.Barrier();
1429  std::vector<std::vector<char> >::iterator rbi = RecvBuffers.begin();
1430  std::vector<int>::iterator rri = recvranks.begin();
1431  // std::vector<std::vector<char> > testbufr(recvranks.size());
1432  while(rbi != RecvBuffers.end()){
1433  // int iiindex = rbi - RecvBuffers.begin();
1434  // testbufr[iiindex].resize(10000);
1435  // PatchComm.ARecv(testbufr[iiindex],*rri);
1436  if(ouf)
1437  *ouf << "Receiving from " << *rri << std::endl;
1438  PatchComm.ARecv(*rbi,*rri);
1439  rbi++;
1440  rri++;
1441  }
1442  if(ouf)
1443  *ouf << "SendBuffers size = " << SendBuffers.size();
1444  std::vector<std::vector<char> >::iterator sbi = SendBuffers.begin();
1445  rri = sendranks.begin();
1446  // std::vector<std::vector<char> > testbufs(sendranks.size());
1447  while(sbi != SendBuffers.end()){
1448  // int iiindex = sbi - SendBuffers.begin();
1449  // testbufs[iiindex].resize(10000);
1450  // PatchComm.ASend(testbufs[iiindex],*rri);
1451  if(ouf)
1452  *ouf << "Sending to " << *rri << std::endl;
1453  PatchComm.ASend(*sbi,*rri);
1454  sbi++;
1455  rri++;
1456  }
1457  PatchComm.WaitAll();
1458  PatchComm.Barrier();
1459  // if(ouf){
1460  // *ouf << "TRAIL: All communication are completed." << std::endl << std::flush;
1461  // }
1462 
1463  if(ouf){
1464  *ouf << "TRAIL: All communication are initiated." << std::endl;
1465  }
1466  // Now all the communcation is initiated. We'll do some
1467  // work while the communication completes.
1468  // Populate the FE representation with the local values
1469  // {
1470  if(true){
1471  // Determine the indices of the FD vertices in the FE
1472  // representation
1473  // std::vector<int> flatindices;
1474  // UnsExtent.GetFlatIndices(FDExtent,flatindices);
1475  PatchComm.Barrier();
1476  if(ouf){
1477  *ouf << "TRAIL: Getting indices of fd nodes in unstructured representation." << std::endl;
1478  std::vector<int> uflat;
1479  std::vector<int> dflat;
1480  UnsExtent.Flatten(uflat);
1481  FDExtent.Flatten(dflat);
1482  *ouf << "TRAIL: finding fd_extent: (["
1483  << dflat[0] << "," << dflat[1] << "],["
1484  << dflat[2] << "," << dflat[3] << "],["
1485  << dflat[4] << "," << dflat[5] << "])"
1486  << std::endl;
1487  *ouf << "TRAIL: in uns_extent: (["
1488  << uflat[0] << "," << uflat[1] << "],["
1489  << uflat[2] << "," << uflat[3] << "],["
1490  << uflat[4] << "," << uflat[5] << "])"
1491  << std::endl;
1492  }
1493  // std::vector<int> flatindices;
1494  // UnsExtent.GetFlatIndices(FDExtent,flatindices);
1495  // if(ouf){
1496  // *ouf << "TRAIL: Made it past finding indices." << std::endl
1497  // << "TRAIL: found overlap: ";
1498  // IRAD::Util::DumpContents(*ouf,theflatindices," ");
1499  // *ouf << std::endl << std::flush;
1500  // }
1501 
1502  PatchComm.Barrier();
1503  if(ouf){
1504  *ouf << "TRAIL: Ready to copy local FD data into FE arrays." << std::endl << std::flush;
1505  }
1506 
1507  // For each attribute, copy the values from the FD vertices
1508  // to the FE vertices
1509  int attindex = 0;
1510  ai = atts.begin();
1511  while(ai != atts.end()){
1512  std::string att = *ai++;
1513  std::string tatt = destwin+"."+att;
1514  std::string satt = srcwin+"."+att;
1515  // First, make sure the destination attribute exists, if not then create and size it
1516  // {
1517  if(COM_get_attribute_handle(tatt) <= 0){
1518  if(ouf && !rank)
1519  *ouf << "TRAIL: " << tatt << " did not exist." << std::endl;
1520  int ncomp = 0;
1521  char loc;
1522  COM_Type type;
1523  std::string unit;
1524  COM_get_attribute(satt,&loc,&type,&ncomp,&unit);
1525  COM_new_attribute(tatt,loc,type,ncomp,unit.c_str());
1526  // COM_set_size(tatt,pane_id,UnsExtent.NNodes());
1527  COM_resize_array(tatt,pane_id);
1528  if(ouf && !rank)
1529  *ouf << "TRAIL: Initializing buffer to 0." << std::endl;
1530  char *trg_ptr = NULL;
1531  COM_get_array(tatt.c_str(),pane_id,&trg_ptr);
1532  if(trg_ptr){
1533  int cnt = 0;
1534  for(int acount = 0;acount < UnsExtent.NNodes();acount++)
1535  for(int ncomp = 0; ncomp < asizes[attindex];ncomp++)
1536  trg_ptr[cnt++] = 0;
1537  }
1538  // COM_window_init_done(destwin,0);
1539  }
1540  else {
1541  if(ouf && !rank)
1542  *ouf << "TRAIL: Attribute appeared to exist on target" << std::endl;
1543  COM_resize_array(tatt,pane_id);
1544  if(ouf && !rank)
1545  *ouf << "TRAIL: Initializing buffer to 0." << std::endl;
1546  char *trg_ptr = NULL;
1547  COM_get_array(tatt.c_str(),pane_id,&trg_ptr);
1548  assert(trg_ptr != NULL);
1549  if(trg_ptr){
1550  int cnt = 0;
1551  for(int acount = 0;acount < UnsExtent.NNodes();acount++)
1552  for(int ncomp = 0; ncomp < asizes[attindex];ncomp++)
1553  trg_ptr[cnt++] = 0;
1554  }
1555  // COM_window_init_done(destwin,0);
1556  }
1557  char *src_ptr = NULL;
1558  char *trg_ptr = NULL;
1559  COM_get_array(tatt.c_str(),pane_id,&trg_ptr);
1560  assert(trg_ptr != NULL);
1561  COM_get_array(satt.c_str(),pane_id,&src_ptr);
1562  assert(src_ptr != NULL);
1563  std::vector<int>::iterator fii = theflatindices.begin();
1564  int count = 0;
1565  while(fii != theflatindices.end()){
1566  int index = *fii++;
1567  memcpy(&trg_ptr[(index-1)*asizes[attindex]],&src_ptr[count++*asizes[attindex]],asizes[attindex]);
1568  }
1569  attindex++;
1570  }
1571  // }
1572  PatchComm.Barrier();
1573  }
1574  // Wait for pending communication
1575  // PatchComm.WaitAll();
1576  // if(ouf && !rank)
1577  // *ouf << "TRAIL: All communication completed." << std::endl;
1578  PatchComm.Barrier();
1579  if(ouf)
1580  *ouf << "TRAIL: Base copy completed." << std::endl << std::flush;
1581  // Populate the FE representation from the received info
1582  if(true){
1583  // Copying out of the receive buffers is dead easy
1584  std::vector<Mesh::BSExtent<int> >::iterator rei = RecvExtents.begin();
1585  int recvindex = 0;
1586  while(rei != RecvExtents.end()){
1587  std::vector<int> flatindices;
1588  unsigned int nnodes = rei->NNodes();
1589  // UnsExtent.GetFlatIndices(*rei++,flatindices);
1590  UnsExtent.GetFlatIndices(*rei,flatindices);
1591  std::vector<int>::iterator fii = flatindices.begin();
1592  assert(nnodes == flatindices.size());
1593  for(unsigned int i = 0;i < nnodes;i++){
1594  int index = *fii++;
1595  int attindex = 0;
1596  ai = atts.begin();
1597  int byte_offset = 0;
1598  while(ai != atts.end()){
1599  std::string tatt = destwin+"."+*ai;
1600  char *trg_ptr = NULL;
1601  COM_get_array(tatt.c_str(),pane_id,&trg_ptr);
1602  assert(trg_ptr != NULL);
1603  // memcpy(&trg_ptr[(index-1)*asizes[attindex]],&RecvBuffers[recvindex][i*bsize+byte_offset],asizes[attindex]);
1604  memcpy(&trg_ptr[(index-1)*asizes[attindex]],
1605  &RecvBuffers[recvindex][i*bsize+byte_offset],asizes[attindex]);
1606  byte_offset += asizes[attindex++];
1607  ai++;
1608  }
1609  }
1610  recvindex++;
1611  rei++;
1612  }
1613  }
1614  } // if(patch_color == 1)
1615  } // if(true) [debugging]
1616  PatchComm.Barrier();
1617  } // loop through patches on this block
1618  } // if(block_color == 1)
1619  } // if(true) [debugging]
1620  BlockComm.Barrier();
1621  } // loop through global blocks
1622  } // if(true) [debugging]
1623  BaseComm.Barrier();
1624  if(ouf && !rank)
1625  *ouf << "TRAIL: All processors done with transfer." << std::endl;
1626  COM_window_init_done(destwin.c_str());
1627  BaseComm.Barrier();
1628  return(0);
1629 }
1630 
1651 int
1652 TRAIL_FD2FE_WinCreate2(const std::string &wname,const std::string &outwname,std::ostream *ouf)
1653 {
1654  // Get the window communicator
1655  MPI_Comm communicator = MPI_COMM_NULL;
1656  COM_get_communicator(wname.c_str(),&communicator);
1657  IRAD::Comm::CommunicatorObject BaseComm(communicator);
1658  int nprocs = BaseComm.Size();
1659  int rank = BaseComm.Rank();
1660  if(ouf)
1661  *ouf << "TRAIL_AddBlockStructuredGhostZone::Nprocs = " << nprocs << std::endl;
1662  std::vector<int> pane_ids;
1663  COM_get_panes(wname.c_str(),pane_ids);
1664  if(ouf)
1665  *ouf << "Number of panes: " << pane_ids.size() << std::endl;
1666  // Form a list of unique block id's across all processors
1667  std::vector<int> local_block_ids;
1668  TRAIL_GetPanelAttribute(wname,"block_id",local_block_ids);
1669  if(ouf){
1670  *ouf << "Local Block_ids: ";
1671  IRAD::Util::DumpContents(*ouf,local_block_ids," ");
1672  *ouf << std::endl;
1673  }
1674  std::vector<int> global_blocks;
1675  TRAIL_UniqueAcrossProcs(local_block_ids,global_blocks,BaseComm.World());
1676  if(ouf){
1677  *ouf << "Global Block_ids: ";
1678  IRAD::Util::DumpContents(*ouf,global_blocks," ");
1679  *ouf << std::endl;
1680  }
1681  // Now all_block_ids is an array of all the unique block ids across all procs
1682  // For each block
1683  // std::vector<int>::iterator bii = global_blocks.begin();
1684  std::string wname2 = outwname; // wname+"_uns";
1685  COM_new_window(wname2.c_str());
1686  COM_new_attribute((wname2+".local_extent").c_str(), 'p',COM_INTEGER,6,"");
1687  COM_new_attribute((wname2+".global_extent").c_str(), 'p',COM_INTEGER,6,"");
1688  COM_new_attribute((wname2+".shared_extents").c_str(),'p',COM_INTEGER,6,"");
1689  COM_new_attribute((wname2+".shared_panes").c_str(), 'p',COM_INTEGER,1,"");
1690  COM_new_attribute((wname2+".bcflag").c_str(), 'p',COM_INTEGER,1,"");
1691  // COM_new_attribute((wname2+".send_extent").c_str(), 'p',COM_INTEGER,6,"");
1692  // COM_new_attribute((wname2+".send_panes").c_str(), 'p',COM_INTEGER,1,"");
1693  // COM_new_attribute((wname2+".recv_panes").c_str(), 'p',COM_INTEGER,1,"");
1694  // COM_new_attribute((wname2+".recv_extent").c_str(), 'p',COM_INTEGER,6,"");
1695  COM_new_attribute((wname2+".block_id").c_str(), 'p',COM_INTEGER,1,"");
1696  COM_new_attribute((wname2+".patch_id").c_str(), 'p',COM_INTEGER,1,"");
1697  // int bcflag = 1;
1698  // COM_set_array((wname2+".bcflag").c_str()),
1699  // COM_new_attribute((wname+".send_extent").c_str(), 'p',COM_INTEGER,6,"");
1700  // COM_new_attribute((wname+".send_panes").c_str(), 'p',COM_INTEGER,1,"");
1701  // COM_new_attribute((wname+".recv_extent").c_str(), 'p',COM_INTEGER,6,"");
1702  // COM_new_attribute((wname+".recv_panes").c_str(), 'p',COM_INTEGER,1,"");
1703 
1704  // Loop through the global block ids on all processors
1705  std::vector<int>::iterator bii = global_blocks.begin();
1706  while(bii != global_blocks.end()){
1707  int block_id = *bii++;
1708  std::vector<int>::iterator fi = std::find(local_block_ids.begin(),local_block_ids.end(),block_id);
1709  int block_color = 0;
1710  if(fi != local_block_ids.end()) // Then this processor has data for the given block
1711  block_color = 1;
1712  // Split the communicator into haves and have nots for this block
1713  if(ouf)
1714  *ouf << "Splitting communicator for blocks." << std::endl;
1715  IRAD::Comm::CommunicatorObject BlockComm;
1716  BaseComm.Split(block_color,rank,BlockComm);
1717  if(block_color == 1){ // all the guys with data on this block
1718  int block_nproc = BlockComm.Size();
1719  int block_rank = BlockComm.Rank();
1720  if(ouf){
1721  *ouf << "Processor " << rank << " has block rank "
1722  << block_rank << "/" << block_nproc << std::endl;
1723  }
1724  // Each block can have multiple "patches". Each patch should be
1725  // contained in its own pane in the source window. confusing naming:
1726  // block = grid, patch = block, new word patch
1727  // Form a list of unique patch id's across all processors
1728  std::vector<int> local_patch_ids;
1729  std::vector<int> panes;
1730  TRAIL_ExtractPanes(wname,"block_id",block_id,panes); // get pane id's of panes with block id = block_id
1731  if(ouf){
1732  *ouf << "Found " << panes.size() << " local panes: ";
1733  IRAD::Util::DumpContents(*ouf,panes," ");
1734  *ouf << std::endl;
1735  }
1736  std::vector<int>::iterator pi = panes.begin();
1737  while(pi != panes.end()){
1738  int *patch_id;
1739  COM_get_array((wname+".patch_id").c_str(),*pi,&patch_id);
1740  if(patch_id)
1741  local_patch_ids.push_back(*patch_id);
1742  pi++;
1743  }
1744  if(ouf){
1745  *ouf << "Found " << local_patch_ids.size() << " local patch_ids: ";
1746  IRAD::Util::DumpContents(*ouf,local_patch_ids," ");
1747  *ouf << std::endl;
1748  }
1749  std::vector<int> all_patches(block_nproc);
1750  TRAIL_UniqueAcrossProcs(local_patch_ids,all_patches,BlockComm.World());
1751  // Now all_patches is an array of all the unique patch id across all procs having the current block
1752  if(ouf){
1753  *ouf << "Found " << all_patches.size() << " global patch_ids: ";
1754  IRAD::Util::DumpContents(*ouf,all_patches," ");
1755  *ouf << std::endl;
1756  }
1757  std::vector<int>::iterator pii = all_patches.begin();
1758  while(pii != all_patches.end()){ // For each patch
1759  int patch_id = *pii++;
1760  // Determine if local processor owns part of the patch
1761  // Split the communicator into haves and have nots
1762  int patch_color = 0;
1763  IRAD::Comm::CommunicatorObject PatchComm;
1764  std::vector<int>::iterator fp = std::find(local_patch_ids.begin(),local_patch_ids.end(),patch_id);
1765  if(fp != local_patch_ids.end())
1766  patch_color = 1;
1767  if(ouf)
1768  *ouf << "Splitting communicator for patches." << std::endl;
1769  BlockComm.Split(patch_color,block_rank,PatchComm);
1770  if(patch_color == 1) { // all of us that have data on the given block/patch
1771  int patch_nproc = PatchComm.Size();
1772  int patch_rank = PatchComm.Rank();
1773  std::vector<int> patch_pane;
1774  TRAIL_ExtractPanes(wname,"patch_id",patch_id,patch_pane); // get the pane for this patch (hopefully only 1)
1775  int *global_patch_extent_ptr = NULL;
1776  COM_get_array((wname+".global_extent").c_str(),patch_pane[0],&global_patch_extent_ptr);
1777  if(!global_patch_extent_ptr){
1778  if(ouf)
1779  *ouf << "ERROR: Window " << wname << " has no global_patch_extent attribute." << std::endl;
1780  std::cerr << "ERROR: Window " << wname << " has no global_patch_extent attribute." << std::endl;
1781  assert(0);
1782  }
1783  if(ouf){
1784  *ouf << "global_patch_extent: (["
1785  << global_patch_extent_ptr[0] << "," << global_patch_extent_ptr[1] << "],["
1786  << global_patch_extent_ptr[2] << "," << global_patch_extent_ptr[3] << "],["
1787  << global_patch_extent_ptr[4] << "," << global_patch_extent_ptr[5] << "])"
1788  << std::endl;
1789  }
1790  int *local_patch_extent_ptr = NULL;
1791  COM_get_array((wname+".local_extent").c_str(),patch_pane[0],&local_patch_extent_ptr);
1792  if(!local_patch_extent_ptr){
1793  if(ouf)
1794  *ouf << "ERROR: Window " << wname << " has no local_patch_extent attribute." << std::endl;
1795  std::cerr << "ERROR: Window " << wname << " has no local_patch_extent attribute." << std::endl;
1796  assert(0);
1797  }
1798  if(ouf){
1799  *ouf << "local_patch_extent: (["
1800  << local_patch_extent_ptr[0] << "," << local_patch_extent_ptr[1] << "],["
1801  << local_patch_extent_ptr[2] << "," << local_patch_extent_ptr[3] << "],["
1802  << local_patch_extent_ptr[4] << "," << local_patch_extent_ptr[5] << "])"
1803  << std::endl;
1804  }
1805  std::vector<int> local_patch_extent;
1806  std::vector<int> flat_local_patch_extents(6*patch_nproc);
1807  IRAD::Util::CopyIntoContainer(local_patch_extent,local_patch_extent_ptr,6);
1808  Mesh::BSExtent<int> LocalPatchExtent(local_patch_extent);
1809  std::vector<int> global_patch_extent;
1810  IRAD::Util::CopyIntoContainer(global_patch_extent,global_patch_extent_ptr,6);
1811  Mesh::BSExtent<int> GlobalPatchExtent(global_patch_extent);
1812 
1813  // Here's a rather simple bit of code that extends the patch to
1814  // to the right if needed.
1815  Mesh::BSExtent<int> ExtendedPatchExtent(LocalPatchExtent);
1816  Mesh::BSExtent<int>::iterator lpei = LocalPatchExtent.begin();
1817  Mesh::BSExtent<int>::iterator gpei = GlobalPatchExtent.begin();
1818  Mesh::BSExtent<int>::iterator epei = ExtendedPatchExtent.begin();
1819  while(lpei != LocalPatchExtent.end()){
1820  if((*lpei)[1] < (*gpei)[1]){
1821  (*epei)[1]++;
1822 
1823  }
1824  lpei++;
1825  gpei++;
1826  epei++;
1827  }
1828 
1829 
1830  // Communicate the extended patch extents to all procs having
1831  // nodes in this patch
1832  std::vector<int> extended_patch_extent;
1833  ExtendedPatchExtent.Flatten(extended_patch_extent);
1834  std::vector<int> all_extended_patch_extents(6*patch_nproc);
1835  // PatchComm.AllGather<std::vector<int>,int>(extended_patch_extent,all_extended_patch_extents,6,6);
1836  PatchComm.AllGather(extended_patch_extent,all_extended_patch_extents,6,6);
1837  // PatchComm.AllGather<std::vector<int>,int>(local_patch_extent,all_local_patch_extents,6,6);
1838  std::vector<int> all_pane_ids(patch_nproc);
1839  // PatchComm.AllGather<std::vector<int>,int>(patch_pane[0],all_pane_ids);
1840  PatchComm.AllGather(patch_pane[0],all_pane_ids);
1841  if(ouf)
1842  *ouf << "All extended patch extents communicated." << std::endl;
1843 
1844  // Fill the search pool, these are the grid extents of all
1845  // remote panes.
1846  std::vector<Mesh::BSExtent<int> > ExtentPool;
1847  for(int pindex = 0;pindex < patch_nproc;pindex++){
1848  if(pindex != patch_rank){
1849  Mesh::BSExtent<int> swimmer(&all_extended_patch_extents[6*pindex]);
1850  ExtentPool.push_back(swimmer);
1851  }
1852  }
1853 
1854  if(ouf)
1855  *ouf << "Search pool filled." << std::endl;
1856  // Find the shared nodes
1857  std::vector<Mesh::BSExtent<int> > shared_extents;
1858  std::vector<int> neighborranks;
1859  std::vector<int> neighborpanes;
1860  if(true){ // Scope temporary vars
1861  std::vector<int> neighborindices;
1862  if(ouf)
1863  *ouf << "Finding shared nodes." << std::endl;
1864  ExtendedPatchExtent.FindSharedNodes(ExtentPool,shared_extents,neighborindices);
1865  if(ouf){
1866  *ouf << "Shared nodes found." << std::endl
1867  << "Neighbor indices: ";
1868  IRAD::Util::DumpContents(*ouf,neighborindices," ");
1869  *ouf << std::endl << "Shared extents: " << std::endl;
1870  std::vector<Mesh::BSExtent<int> >::iterator bsmi = shared_extents.begin();
1871  while(bsmi != shared_extents.end()){
1872  unsigned int nid = bsmi - shared_extents.begin();
1873  *ouf << "Shared extent " << nid << ":" << std::endl;
1874  std::vector<int> tempflat_extent;
1875  bsmi++->Flatten(tempflat_extent);
1876  *ouf << "shared_extent: (["
1877  << tempflat_extent[0] << "," << tempflat_extent[1] << "],["
1878  << tempflat_extent[2] << "," << tempflat_extent[3] << "],["
1879  << tempflat_extent[4] << "," << tempflat_extent[5] << "])"
1880  << std::endl;
1881  }
1882  }
1883  std::vector<int>::iterator nii = neighborindices.begin();
1884  while(nii != neighborindices.end()){
1885  int nrank = *nii++;
1886  if(nrank >= patch_rank)
1887  nrank++;
1888  neighborranks.push_back(nrank);
1889  neighborpanes.push_back(all_pane_ids[nrank]);
1890  }
1891  }
1892  if(ouf){
1893  *ouf << "Neighbor ranks: ";
1894  IRAD::Util::DumpContents(*ouf,neighborranks," ");
1895  *ouf << std::endl << "Neighborpanes: ";
1896  IRAD::Util::DumpContents(*ouf,neighborpanes," ");
1897  *ouf << std::endl;
1898  }
1899 
1900  if(ouf)
1901  *ouf << "Copying attributes." << std::endl;
1902  TRAIL_Copy2Attribute((wname2+".local_extent"),extended_patch_extent,patch_pane[0]);
1903  TRAIL_Copy2Attribute((wname2+".global_extent"),global_patch_extent,patch_pane[0]);
1904  TRAIL_SetAttribute((wname2+".block_id"),patch_pane[0],block_id);
1905  TRAIL_SetAttribute((wname2+".patch_id"),patch_pane[0],patch_id);
1906  Mesh::Connectivity conn;
1907  ExtendedPatchExtent.Sync();
1908  if(ouf)
1909  *ouf << "Creating unstructured mesh." << std::endl;
1910  ExtendedPatchExtent.CreateUnstructuredMesh(conn);
1911  conn.Sync();
1912  conn.SyncSizes();
1913  bool flip = false;
1914  int flip_hndl = COM_get_attribute_handle((wname+".flip").c_str());
1915  if(flip_hndl > 0){
1916  int *flip_ptr = NULL;
1917  COM_get_array((wname+".flip").c_str(),patch_pane[0],&flip_ptr);
1918  if(flip_ptr != NULL)
1919  flip = (*flip_ptr > 0);
1920  }
1921  if(flip){
1922  if(ouf)
1923  *ouf << "Flipping unstructured surface." << std::endl;
1924  Mesh::Connectivity::iterator connit = conn.begin();
1925  while(connit != conn.end())
1926  {
1927  Mesh::GenericCell_2 ge(connit->size());
1928  ge.ReOrient(*connit++);
1929  }
1930  }
1931 // if(ouf){
1932 // std::vector<int> tempout;
1933 // ExtendedPatchExtent.Flatten(tempout);
1934 // *ouf << "Unstructured extent: ";
1935 // IRAD::Util::DumpContents(*ouf,tempout," ");
1936 // *ouf << std::endl;
1937 // *ouf << "Connectivity: " << std::endl
1938 // << conn << std::endl;
1939 // }
1940  if(ouf)
1941  *ouf << "Setting window data." << std::endl;
1942  std::vector<int> flatcon;
1943  unsigned int nnodes_new = ExtendedPatchExtent.NNodes();
1944  conn.Flatten(flatcon);
1945  COM_set_size((wname2+".nc").c_str(),patch_pane[0],nnodes_new);
1946  COM_resize_array((wname2+".nc").c_str(),patch_pane[0]);
1947  COM_set_size((wname2+".:q4:").c_str(),patch_pane[0],conn.Nelem());
1948  COM_resize_array((wname2+".:q4:").c_str(),patch_pane[0]);
1949  COM_set_size((wname2+".bcflag").c_str(),patch_pane[0],1);
1950  COM_resize_array((wname2+".bcflag").c_str(),patch_pane[0]);
1951  int *bcflag;
1952  COM_get_array((wname2+".bcflag").c_str(),patch_pane[0],&bcflag);
1953  *bcflag = 1;
1954 
1955  int *conndat;
1956  COM_get_array((wname2+".:q4:").c_str(),patch_pane[0],&conndat);
1957  memcpy(conndat,&flatcon[0],sizeof(int)*4*conn.Nelem());
1958  std::vector<std::vector<int> > flattened_shared_extents;
1959  std::vector<Mesh::BSExtent<int> >::iterator sei = shared_extents.begin();
1960  while(sei != shared_extents.end()){
1961  std::vector<int> flatextent;
1962  sei++->Flatten(flatextent);
1963  flattened_shared_extents.push_back(flatextent);
1964  }
1965  TRAIL_Copy2Attribute((wname2+".shared_panes"),neighborpanes,patch_pane[0],neighborpanes.size());
1966  TRAIL_Copy2Attribute((wname2+".shared_extents"),flattened_shared_extents,patch_pane[0]);
1967  if(ouf)
1968  *ouf << "Building pconn." << std::endl;
1969  std::vector<int> pconn;
1970  pconn.push_back(neighborpanes.size());
1971  std::vector<int> sortedpanes(neighborpanes);
1972  std::sort(sortedpanes.begin(),sortedpanes.end());
1973  std::vector<int>::iterator spit = sortedpanes.begin();
1974  while(spit != sortedpanes.end()){
1975  int remote_pane_id = *spit++;
1976  pconn.push_back(remote_pane_id);
1977  int remote_pane_index = std::find(neighborpanes.begin(),neighborpanes.end(),remote_pane_id) - neighborpanes.begin();
1978  std::vector<int> sharednodes;
1979  ExtendedPatchExtent.GetFlatIndices(shared_extents[remote_pane_index],sharednodes);
1980  pconn.push_back(sharednodes.size());
1981  std::vector<int>::iterator sni = sharednodes.begin();
1982  while(sni != sharednodes.end())
1983  pconn.push_back(*sni++);
1984  }
1985  TRAIL_Copy2Attribute((wname2+".pconn"),pconn,patch_pane[0],pconn.size());
1986  if(ouf){
1987  *ouf << "PConn: ";
1988  IRAD::Util::DumpContents(*ouf,pconn," ");
1989  *ouf << std::endl;
1990  }
1991  }
1992  PatchComm.Barrier();
1993  }
1994  }
1995  BlockComm.Barrier();
1996  }
1997  if(ouf)
1998  *ouf << "window done, roccom finalizing it." << std::endl;
1999  COM_window_init_done(wname2.c_str());
2000  if(ouf)
2001  *ouf << "window finalized." << std::endl;
2002  if(ouf)
2003  *ouf << "end barrier" << std::endl;
2004  BaseComm.Barrier();
2005  if(ouf)
2006  *ouf << "all processors exiting." << std::endl;
2007  return(0);
2008 }
2009 
2010 void TRAIL_Window2UnstructuredMesh(const std::string &wname,std::vector<Mesh::UnstructuredMesh> &meshes,
2011  std::vector<SolnMetaData> &smdv,
2012  std::vector<std::vector<std::vector<double> > > &soln_data,
2013  int verblevel, bool no_ghost)
2014 {
2015  std::vector<int> pane_ids;
2016  std::vector<std::vector<std::vector<double> > > temp_soln;
2017  TRAIL_GetWindowSolnMetaData(wname,smdv,verblevel>1);
2018  TRAIL_GetWindowSolnData(wname,temp_soln,smdv,verblevel>1);
2019  soln_data.resize(temp_soln.size());
2020  meshes.resize(temp_soln.size());
2021  COM_get_panes(wname.c_str(),pane_ids);
2022  if(verblevel > 1)
2023  std::cout << "TRAIL_Window2UnstructuredMesh::Processing " << pane_ids.size()
2024  << " panes." << std::endl;
2025  std::vector<int>::iterator pii = pane_ids.begin();
2026  while(pii != pane_ids.end()){
2027  int pane_index = pii - pane_ids.begin();
2028  int pane_id = *pii++;
2029  const int *conn_ptr = NULL;
2030  double *coordinates_ptr = NULL;
2031  COM_get_array_const((wname+".nc").c_str(),pane_id,(const double **)&coordinates_ptr);
2032  int nreal = 0;
2033  int nghost = 0;
2034  COM_get_array_const((wname+".:st3:").c_str(),pane_id,&conn_ptr);
2035  int nnodes = 0;
2036  if(conn_ptr){
2037  COM_get_size((wname+".:st3:").c_str(),pane_id,&nreal,&nghost);
2038  if(verblevel)
2039  std::cout << "Block Structured pane " << pane_id << ": " << std::endl
2040  << "Coordinate array size (i,j,k): (" << conn_ptr[0]
2041  << "," << conn_ptr[1] << "," << conn_ptr[2] << ")" << std::endl
2042  << "Ghost zone is " << nghost << " cells wide."
2043  << std::endl;
2044 
2045 
2046  nnodes = conn_ptr[0]*conn_ptr[1]*conn_ptr[2];
2047  meshes[pane_index].nc.init_copy(nnodes,coordinates_ptr);
2048  std::vector<Mesh::IndexType> extent;
2049  for(int i = 0;i < 3;i++){
2050  extent.push_back(1);
2051  extent.push_back(conn_ptr[i]);
2052  }
2053  // resize soln_data to the number of attributes
2054  soln_data[pane_index].resize(temp_soln[pane_index].size());
2055 
2056  Mesh::BSExtent<Mesh::IndexType> bsextent(extent);
2057  if(no_ghost){
2058  if(verblevel > 1)
2059  std::cout << "TRAIL_Window2UnstructuredMesh::Stripping ghosts from block structured Window"
2060  << std::endl;
2061  std::vector<Mesh::IndexType> real_extent;
2062  for(int i = 0;i < 3;i++){
2063  real_extent.push_back(1+nghost);
2064  real_extent.push_back(conn_ptr[i]-nghost);
2065  }
2066  if(verblevel > 1)
2067  std::cout << "TRAIL_Window2UnstructuredMesh::Getting real node indices"
2068  << std::endl;
2069  Mesh::BSExtent<Mesh::IndexType> realextent(real_extent);
2070  std::vector<Mesh::IndexType> real_indices;
2071  bsextent.GetFlatIndices(realextent,real_indices);
2072  unsigned int nreal_nodes = real_indices.size();
2073  std::vector<double> real_coordinates(nreal_nodes*3);
2074  std::vector<Mesh::IndexType>::iterator rii = real_indices.begin();
2075  Mesh::IndexType rni = 0;
2076  if(verblevel > 1)
2077  std::cout << "TRAIL_Window2UnstructuredMesh::Extracting real coordinates"
2078  << std::endl;
2079  while(rii != real_indices.end()){
2080  Mesh::IndexType cindex = *rii++;
2081  double *crdptr = meshes[pane_index].nc[cindex];
2082  real_coordinates[rni++] = *crdptr++;
2083  real_coordinates[rni++] = *crdptr++;
2084  real_coordinates[rni++] = *crdptr;
2085  }
2086  if(verblevel > 1)
2087  std::cout << "TRAIL_Window2UnstructuredMesh::Replacing coordinates with real"
2088  << std::endl;
2089  meshes[pane_index].nc.init_copy(nreal_nodes,&real_coordinates[0]);
2090  if(verblevel > 1)
2091  std::cout << "TRAIL_Window2UnstructuredMesh::Real nodal coordinates extracted."
2092  << std::endl;
2093  for(int i = 0;i < 3;i++){
2094  real_extent[i*2] = 1;
2095  real_extent[i*2+1] -= nghost;
2096  }
2097  Mesh::BSExtent<Mesh::IndexType> rext(real_extent);
2098  if(verblevel > 1)
2099  std::cout << "TRAIL_Window2UnstructuredMesh::Creating mesh with real vertices/cells"
2100  << std::endl;
2101  rext.CreateUnstructuredMesh(meshes[pane_index].con);
2102  meshes[pane_index].con.Sync();
2103  meshes[pane_index].con.SyncSizes();
2104  if(verblevel > 1)
2105  std::cout << "TRAIL_Window2UnstructuredMesh::Nodes: " << meshes[pane_index].nc.size()
2106  << ", Cells: " << meshes[pane_index].con.Nelem() << std::endl
2107  << "TRAIL_Window2UnstructuredMesh::Extracting real cell indices."
2108  << std::endl;
2109 
2110  std::vector<Mesh::IndexType> cell_extent;
2111  for(int i = 0;i < 3;i++){
2112  cell_extent.push_back(1);
2113  cell_extent.push_back(conn_ptr[i]-1);
2114  }
2115  Mesh::BSExtent<Mesh::IndexType> cellextent(cell_extent);
2116  std::vector<Mesh::IndexType> real_cell_extent;
2117  for(int i = 0;i < 3;i++){
2118  real_cell_extent.push_back(1+nghost);
2119  real_cell_extent.push_back(conn_ptr[i]-1-nghost);
2120  }
2121  Mesh::BSExtent<Mesh::IndexType> realcellextent(real_cell_extent);
2122  std::vector<Mesh::IndexType> real_cell_indices;
2123  cellextent.GetFlatIndices(realcellextent,real_cell_indices);
2124 
2125  if(verblevel > 1)
2126  std::cout << "TRAIL_Window2UnstructuredMesh::Extracting real solution data"
2127  << std::endl;
2128  std::vector<std::vector<double> >::iterator tsai = temp_soln[pane_index].begin();
2129  std::vector<std::vector<double> >::iterator sdai = soln_data[pane_index].begin();
2130  while(tsai != temp_soln[pane_index].end()){
2131  unsigned int att_index = tsai - temp_soln[pane_index].begin();
2132  unsigned int ncomp = smdv[att_index].ncomp;
2133  if(!tsai->empty()){ // make sure there's data for this attribute
2134  if(smdv[att_index].loc == 'N' || smdv[att_index].loc == 'n'){
2135  // It's a nodal attribute
2136  sdai->resize(real_indices.size()*ncomp);
2137  if(verblevel > 1)
2138  std::cout << "TRAIL_Window2UnstructuredMesh::Processing "
2139  << ncomp << " nodal components of "
2140  << smdv[att_index].name << std::endl;
2141  std::vector<double>::iterator sddi = sdai->begin();
2142  rii = real_indices.begin();
2143  while(rii != real_indices.end()){
2144  Mesh::IndexType real_index = *rii++ - 1;
2145  real_index *= ncomp;
2146  for(unsigned int ncc = 0;ncc < ncomp;ncc++)
2147  *sddi++ = temp_soln[pane_index][att_index][real_index++];
2148 
2149  }
2150  }
2151  else if(smdv[att_index].loc == 'E' || smdv[att_index].loc == 'e'){
2152  // It's a cellular attribute
2153  sdai->resize(real_cell_indices.size()*ncomp);
2154  std::vector<double>::iterator sddi = sdai->begin();
2155  if(verblevel > 1)
2156  std::cout << "TRAIL_Window2UnstructuredMesh::Processing "
2157  << ncomp << " cellular components of "
2158  << smdv[att_index].name << std::endl;
2159  rii = real_cell_indices.begin();
2160  while(rii != real_cell_indices.end()){
2161  Mesh::IndexType real_index = *rii++ - 1;
2162  real_index *= ncomp;
2163  for(unsigned int ncc = 0;ncc < ncomp;ncc++)
2164  *sddi++ = temp_soln[pane_index][att_index][real_index++];
2165 
2166  }
2167  }
2168  }
2169  else
2170  sdai->resize(0);
2171  tsai->resize(0);
2172  sdai++;
2173  tsai++;
2174  }
2175  temp_soln[pane_index].resize(0);
2176  }
2177  else{
2178  bsextent.CreateUnstructuredMesh(meshes[pane_index].con);
2179  std::vector<std::vector<double> >::iterator tsai = temp_soln[pane_index].begin();
2180  std::vector<std::vector<double> >::iterator sdai = soln_data[pane_index].begin();
2181  while(tsai != temp_soln[pane_index].end()){
2182  // unsigned int att_index = tsai - temp_soln[pane_index].begin();
2183  // unsigned int ncomp = smdv[att_index].ncomp;
2184  sdai->resize(tsai->size());
2185  std::vector<double>::iterator sddi = sdai->begin();
2186  std::vector<double>::iterator tddi = tsai->begin();
2187  while(tddi != tsai->end())
2188  *sddi++ = *tddi++;
2189  tsai->resize(0);
2190  tsai++;
2191  sdai++;
2192  }
2193  temp_soln[pane_index].resize(0);
2194  }
2195  // meshes.push_back(mesh);
2196  }
2197  else{
2198  if(verblevel > 1)
2199  std::cout << "TRAIL_Window2UnstructuredMesh::Processing unstructured window" << std::endl;
2200  // its unstructured - easier to deal with (i.e. it's already
2201  // unstructured!
2202  }
2203  // COM_get_array_cont((
2204  // << "ST3 for geopane " << pane_id << std::endl
2205  // << geo_array_ptr[0] << "," << geo_array_ptr[1]
2206  // << "," << geo_array_ptr[2] << std::endl;
2207  }
2208 // std::vector<Mesh::UnstructuredMesh>::iterator mi = meshes.begin();
2209 // std::ofstream Ouf;
2210 // Ouf.open("test_coords");
2211 // while(mi != meshes.end()){
2212 // GeoPrim::CBox mesh_box;
2213 // GeoPrim::CBox small_box;
2214 // GeoPrim::CBox large_box;
2215 // Ouf << mi->nc << std::endl;
2216 // Mesh::GetMeshBoxes(mi->nc,mi->con,mesh_box,small_box,large_box);
2217 // std::cout << "Mesh BOXES:" << std::endl
2218 // << "Mesh: " << mesh_box << std::endl
2219 // << "Small: " << small_box << std::endl
2220 // << "Large: " << large_box << std::endl;
2221 // mi++;
2222 // }
2223 // Ouf.close();
2224 }
2225 //#endif
2226 
2227 // Serial Function:
2228 // Read the HDF file fname into the window wname
2229 void TRAIL_HDF2Window( const std::string &fname, const std::string &wname,int verb) {
2230  COM_new_window( wname.c_str());
2231  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocin, "IN");
2232  int IN_read;
2233  IN_read = COM_get_function_handle( "IN.read_window");
2234  MPI_Comm comm_null = MPI_COMM_NULL;
2235  std::string bufwin("bufwin");
2236  if(verb)
2237  std::cout << "Reading file " << fname << "..." << std::endl;
2238  COM_call_function( IN_read, fname.c_str(), bufwin.c_str(), &comm_null);
2239  if(verb)
2240  std::cout << "Done reading file " << fname << "." << std::endl;
2241  int IN_obtain = COM_get_function_handle( "IN.obtain_attribute");
2242  int buf_all = COM_get_attribute_handle((bufwin+".all").c_str());
2243  COM_call_function( IN_obtain, &buf_all, &buf_all);
2244  COM_UNLOAD_MODULE_STATIC_DYNAMIC( Rocin, "IN");
2245  if(verb)
2246  std::cout << "Obtained temp window from file " << fname
2247  << ", cloning.." << std::endl;
2248  COM_window_init_done(bufwin.c_str());
2249  // Change the memory layout to contiguous
2250  COM_clone_attribute( (wname+".all").c_str(), (bufwin+".all").c_str(), 1);
2251  COM_window_init_done(wname.c_str());
2252  COM_delete_attribute((bufwin+".atts").c_str());
2253  COM_delete_window( bufwin.c_str());
2254  if(verb)
2255  std::cout << "Window " << wname << " created." << std::endl;
2256 
2257 }
2258 
2259 void
2260 TRAIL_GetWindowSolnMetaData(const std::string &wname,std::vector<SolnMetaData> &smdv,int verblevel)
2261 {
2262  int na = 0;
2263  std::string atts;
2264  COM_get_attributes(wname.c_str(),&na,atts);
2265  int count = 0;
2266  if(na > 0){
2267  std::istringstream Istr(atts);
2268  std::string aname;
2269  while(Istr >> aname) {
2270  std::string rstring(wname+"."+aname);
2271  char loc;
2272  int ncomp;
2273  COM_Type ta;
2274  std::string unit;
2275  COM_get_attribute(rstring.c_str(),&loc,&ta,&ncomp,&unit);
2276  if( (loc == 'E' || loc == 'e' || loc == 'n' || loc == 'N') && ta == COM_DOUBLE)
2277  count++;
2278  }
2279  }
2280  smdv.resize(count);
2281  count = 0;
2282  if(na > 0) {
2283  std::istringstream Istr(atts);
2284  std::string aname;
2285  while(Istr >> aname){
2286  std::string rstring(wname+"."+aname);
2287  char loc;
2288  int ncomp;
2289  COM_Type ta;
2290  std::string unit;
2291  COM_get_attribute(rstring.c_str(),&loc,&ta,&ncomp,&unit);
2292  if( (loc == 'E' || loc == 'e' || loc == 'n' || loc == 'N') && ta == COM_DOUBLE) {
2293  // SolnMetaData smd;
2294  smdv[count].loc = loc;
2295  smdv[count].name = aname;
2296  smdv[count].unit = unit;
2297  smdv[count].ncomp = ncomp;
2298  // smdv[count].push_back(smd);
2299  count++;
2300  }
2301  }
2302  }
2303 }
2304 
2305 void TRAIL_GetWindowSolnData(const std::string &wname,
2306  std::vector<std::vector<std::vector<double> > > &soln_data,
2307  std::vector<SolnMetaData> &smdv,int verblevel)
2308 {
2309 
2310 
2311  int nAttrs;
2312  char *attrStr;
2313  COM_get_attributes(wname.c_str(), &nAttrs, &attrStr);
2314  std::cout << "Number of attributes: " << nAttrs << std::endl;
2315  std::cout << "Attribute string: " << attrStr << std::endl;
2316 
2317  std::vector<int> pane_ids;
2318  COM_get_panes(wname.c_str(),pane_ids);
2319 
2320  if( soln_data.size( ) != pane_ids.size( ) )
2321  soln_data.resize(pane_ids.size() );
2322 
2323  std::vector<int>::iterator pii = pane_ids.begin();
2324  while(pii != pane_ids.end()) {
2325 
2326  int pane_index = pii - pane_ids.begin();
2327  int pane_id = *pii++;
2328 
2329  if(verblevel > 0) {
2330  std::cout << "TRAIL_GetWindowSolnData::Pane id: " << pane_id << std::endl
2331  << "TRAIL_GetWindowSolnData::Number of atts: "
2332  << smdv.size() << std::endl;
2333  }
2334 
2335  soln_data[pane_index].resize(smdv.size());
2336 
2337  //BugFix: update smdi iterator for each pane (gzagaris)
2338  std::vector<SolnMetaData>::iterator smdi = smdv.begin();
2339  while(smdi != smdv.end()) {
2340 
2341  std::string name = wname + "." + smdi ->name;
2342 // std::string name = wname + ".siVel";
2343  unsigned int aindex = smdi - smdv.begin();
2344  int array_size = 0;
2345 
2346  COM_get_size( name.c_str(), pane_id, &array_size);
2347  if(verblevel)
2348  std::cout << smdi->name << " size = " << array_size << std::endl;
2349 
2350  const double *solver_array = NULL;
2351 
2352  COM_get_array_const( name.c_str( ), pane_id, &solver_array);
2353 
2354 // if( array_size > 0 )
2355 // {
2356 // assert( solver_array != NULL );
2357 // }
2358 
2359  if( solver_array ) {
2360 
2361  if(verblevel)
2362  std::cout << "TRAIL_GetWindowSolnData::Copying " << array_size*smdi->ncomp << " doubles into " << smdi->name << " array." << std::endl;
2363 
2364  soln_data[pane_index][aindex].resize(array_size*smdi->ncomp);
2365  std::memcpy(&soln_data[pane_index][aindex][0],solver_array,sizeof(double)*array_size*smdi->ncomp);
2366  std::vector<double>::iterator ai = soln_data[pane_index][aindex].begin();
2367  }
2368  else {
2369 
2370  if(verblevel)
2371  std::cout << "TRAIL_GetWindowSolnData::No " << smdi->name << " array found." << std::endl;
2372 
2373  soln_data[pane_index][aindex].resize(0);
2374  }
2375 
2376  smdi++;
2377 
2378  } // while smdi != smdv.end()
2379 
2380  } // For all panes
2381 
2382 }
2383 
2384 // Read the Rocin control file fname into he window wname on
2385 // processors in comm (or MPI_COMM_NULL for 1 proc) - read only
2386 // panes with matching bcflags (or all panes if bcflag.empty(),
2387 // and if(apply_disp) then apply the uhat displacement field,
2388 // get all attributes if(all) or only the mesh if(!all),
2389 // and read ghost info or not.
2390 void
2391 TRAIL_File2Window( const std::string &fname,
2392  const std::string &wname,
2393  std::vector<int> &bcflags,
2394  MPI_Comm comm,
2395  bool apply_disp,
2396  bool all,
2397  bool with_ghost)
2398 {
2399 
2400  COM_new_window( wname.c_str());
2401  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocin, "TRAILIN");
2402  int IN_read;
2403  IN_read = COM_get_function_handle( "TRAILIN.read_by_control_file");
2404 
2405 
2406  std::string bufwin("bufwin");
2407  COM_call_function( IN_read, fname.c_str(), bufwin.c_str(), &comm);
2408  int IN_obtain = COM_get_function_handle( "TRAILIN.obtain_attribute");
2409 
2410  // If bcflags specified, keep only panes with matching bcflags
2411  if(!bcflags.empty()){
2412  int bcflag = COM_get_attribute_handle((bufwin+".bcflag").c_str());
2413  if (bcflag > 0) {
2414  COM_call_function( IN_obtain, &bcflag, &bcflag);
2415  int npanes, *pane_ids;
2416  COM_get_panes( bufwin.c_str(), &npanes, &pane_ids);
2417 
2418  // remove panes with bcflag not found in bcflags
2419  for ( int i=0; i<npanes; ++i) {
2420  int *flag;
2421  COM_get_array( (bufwin+".bcflag").c_str(), pane_ids[i], &flag);
2422  if ( flag==NULL )
2423  COM_delete_pane( bufwin.c_str(),pane_ids[i]);
2424  bool delite = true;
2425  std::vector<int>::iterator bcfi = bcflags.begin();
2426  while(bcfi != bcflags.end() && delite)
2427  if(*bcfi++ == *flag)
2428  delite = false;
2429  if(delite)
2430  COM_delete_pane( bufwin.c_str(), pane_ids[i]);
2431  }
2432  COM_free_buffer( &pane_ids);
2433  }
2434  }
2435  if(apply_disp){
2436  // This is NOT correct for problems with regression.
2437  int disp_hndl = COM_get_attribute_handle((bufwin+".uhat").c_str());
2438  if(disp_hndl > 0){
2439  COM_call_function( IN_obtain, &disp_hndl, &disp_hndl);
2440  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocblas, "BLAS");
2441  int add = COM_get_function_handle( "BLAS.add");
2442  COM_call_function(IN_obtain,&disp_hndl,&disp_hndl);
2443  int nc_hndl = COM_get_attribute_handle( bufwin + ".nc");
2444  COM_call_function( add, &disp_hndl, &nc_hndl, &nc_hndl);
2445  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocblas,"BLAS");
2446  }
2447  }
2448  // Read in the mesh.
2449  int buf_atts;
2450  if(all)
2451  buf_atts = COM_get_attribute_handle((bufwin+".all").c_str());
2452  else
2453  buf_atts = COM_get_attribute_handle((bufwin+".mesh").c_str());
2454  COM_call_function( IN_obtain, &buf_atts, &buf_atts);
2455  COM_UNLOAD_MODULE_STATIC_DYNAMIC( Rocin, "TRAILIN");
2456 
2457  if(!all)
2458  // Remove all attributes except for the mesh
2459  COM_delete_attribute( (bufwin+".atts").c_str());
2460 
2461  // Change the memory layout to contiguous.
2462  if(all)
2463  COM_clone_attribute( (wname+".all").c_str(),
2464  (bufwin+".all").c_str(), (with_ghost ? 1 : 0));
2465  else
2466  COM_clone_attribute( (wname+".mesh").c_str(),
2467  (bufwin+".mesh").c_str(), (with_ghost ? 1 : 0));
2468 
2469  COM_delete_window( bufwin.c_str());
2470  COM_window_init_done(wname);
2471 }
2472 
2473 
2474 // DO NOT CALL IN PARALLEL: This function processes files for every processor
2475 void
2476 TRAIL_MergeRocinFiles(const std::string &srcname,
2477  const std::string &trgname,
2478  const std::string &path,
2479  double t,unsigned int np,
2480  std::ostream *ouf)
2481 {
2482  std::string homedir(TRAIL_CWD());
2483  std::string timestring(TRAIL_TimeString(t));
2484  std::string filepre(srcname+"_"+timestring+"_");
2485  std::ostringstream Ofstr;
2486  unsigned int id = 0;
2487  if(ouf)
2488  *ouf << "TRAIL_MergeRocinFiles: Entry" << std::endl;
2489  TRAIL_CD(path,ouf);
2490  if(ouf)
2491  *ouf << "Searching for files with prefix: " << filepre << std::endl;
2492  while(id <= np){
2493  std::ifstream Inf;
2494  std::ostringstream Ostr;
2495  Ostr << filepre << std::setw(5) << std::setfill('0')
2496  << id++;
2497  std::string filename(Ostr.str()+"_in.txt");
2498  Inf.open(filename.c_str());
2499  if(Inf){
2500  Ofstr << Inf.rdbuf() << std::endl;
2501  Inf.close();
2502  unlink(filename.c_str());
2503  }
2504  }
2505  std::ofstream Ouf;
2506  Ouf.open((trgname+"_in_"+timestring+".txt").c_str());
2507  Ouf << Ofstr.str();
2508  Ouf.close();
2509  TRAIL_CD(homedir,ouf);
2510  if(ouf)
2511  *ouf << "TRAIL_MergeRocinFiles: Exit" << std::endl;
2512 }
2513 
2514 void
2515 TRAIL_CreateRobustFC(const std::string &wname,const std::string &path)
2516 {
2517  std::ofstream Ouf;
2518  Ouf.open((path+"/"+wname+".fc").c_str());
2519  Ouf << "0.5 0.6 3 0.17365" << std::endl
2520  << "1.3962634 0.314159265 3" << std::endl
2521  << "0.0 0.5 3" << std::endl
2522  << "6 1 1 0" << std::endl
2523  << "2" << std::endl;
2524  Ouf.close();
2525 }
2526 void
2527 TRAIL_CreateRobustFC_old(const std::string &wname,const std::string &path)
2528 {
2529  std::ofstream Ouf;
2530  Ouf.open((path+"/"+wname+".fc").c_str());
2531  Ouf << "0.76604444 0.98480775 3 0.98480775" << std::endl
2532  << "1.3962634 0.314159265 3" << std::endl
2533  << "0.17364818 0.96592583 3" << std::endl
2534  << "6 1 1 0" << std::endl
2535  << "2" << std::endl;
2536  Ouf.close();
2537 }
2538 
2539 // serial function: **do not call on multiple procs**
2540 // Reads inputs for the source and target windows, whos names are
2541 // spec'd by src,trg from the srcpath and trgpath, and writes the common
2542 // refinement to the destpath.
2543 void
2544 TRAIL_AutoSurfer(const std::string &src, const std::string &trg,
2545  const std::string &srcpath, const std::string &trgpath,
2546  const std::string &destpath, double t, MPI_Comm comm,
2547  std::ostream *ouf)
2548 {
2549  std::string timestring(TRAIL_TimeString(t));
2550  std::string srcfile(src + "_in_" + timestring + ".txt");
2551  std::string trgfile(trg + "_in_" + timestring + ".txt");
2552  std::string trailwin(src+"_trail");
2553  std::string homedir(TRAIL_CWD());
2554  std::string format("HDF");
2555  std::string newpath;
2556  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocface, "RFC");
2557  int RFC_overlay = COM_get_function_handle( "RFC.overlay");
2558  int RFC_write = COM_get_function_handle( "RFC.write_overlay");
2559  TRAIL_CreateRobustFC(trg,destpath);
2560  TRAIL_CreateRobustFC(trailwin,destpath);
2561  COM_set_default_communicator(MPI_COMM_NULL);
2562  std::vector<int> bcflags(3);
2563  bcflags[0] = 0;
2564  bcflags[1] = 1;
2565  bcflags[2] = 2;
2566  TRAIL_CD(srcpath,ouf);
2567  if(ouf)
2568  *ouf << "TRAIL_AutoSurfer:: homedir = " << homedir << std::endl
2569  << "TRAIL_AutoSurfer:: CWD = " << TRAIL_CWD() << std::endl
2570  << "TRAIL_AutoSurfer:: Creating common refinement." << std::endl
2571  << "TRAIL_AutoSurfer:: Reading in source surface from " << srcpath
2572  << std::endl;
2573  if(ouf)
2574  std::cout << "Roctrail> Reading in source surface from "
2575  << srcpath << "(test of patience)" << std::endl;
2576  TRAIL_File2Window(srcfile,trailwin,bcflags,MPI_COMM_NULL,false,false,false);
2577  newpath.assign(homedir+"/"+trgpath);
2578  TRAIL_CD(newpath,ouf);
2579  if(ouf)
2580  *ouf << "TRAIL_AutoSurfer: Reading target surface from " << newpath
2581  << std::endl;
2582  if(ouf)
2583  std::cout << "Roctrail> Reading target surface from " << newpath
2584  << std::endl;
2585  TRAIL_File2Window(srcfile,trg,bcflags,MPI_COMM_NULL,false,false,false);
2586  newpath.assign(homedir+"/"+destpath);
2587  TRAIL_CD(newpath,ouf);
2588  int src_mesh = COM_get_attribute_handle( (trailwin+".mesh").c_str());
2589  int trg_mesh = COM_get_attribute_handle( (trg+".mesh").c_str());
2590  if(ouf){
2591  *ouf << "TRAIL_AutoSurfer: Calling Rocface Overlay...."
2592  << std::endl;
2593  std::cout << "Roctrail> Rocface performing overlay..." << std::endl;
2594  }
2595  COM_call_function( RFC_overlay, &src_mesh, &trg_mesh);
2596  if(ouf){
2597  std::cout << "Roctrail> Overlay complete, writing overlay to " << newpath
2598  << std::endl;
2599  *ouf << "TRAIL_AutoSurfer: Writing overlay to " << newpath << "..."
2600  << std::endl;
2601  }
2602  COM_call_function( RFC_write, &src_mesh, &trg_mesh,
2603  trailwin.c_str(), trg.c_str(), format.c_str());
2604  TRAIL_CD(homedir,ouf);
2605  if(ouf)
2606  *ouf << "TRAIL_AutoSurfer: Done. CWD = " << TRAIL_CWD() << std::endl;
2607  COM_delete_window(trg);
2608  COM_delete_window(trailwin);
2609  COM_set_default_communicator(comm);
2610  COM_UNLOAD_MODULE_STATIC_DYNAMIC( Rocface, "RFC");
2611 }
2612 
2613 void TRAIL_FixRocstarData(const std::string &wname,std::ostream *ouf = NULL);
2614 void TRAIL_ExtractSurf0(const std::string &srcwin,
2615  const std::string &trgwin,
2616  std::ostream *ouf = NULL);
2617 
2618 bool
2619 TRAIL_TransferSurfDataFILE
2620 (
2621  const std::string &src,
2622  const std::string &trg,
2623  const std::string &dest,
2624  const std::string &srcpath,
2625  const std::string &trgpath,
2626  const std::string &destpath,
2627  const std::string &crpath,
2628  double t,unsigned int id,
2629  MPI_Comm comm,
2630  std::ostream *ouf)
2631 {
2632  std::string timestring(TRAIL_TimeString(t));
2633  std::string suffix("_in_"+timestring+".txt");
2634  std::string srcfile(src + suffix);
2635  std::string trgfile(trg + suffix);
2636  std::string r_trg(dest);
2637  std::string trailwin(src+"_trail");
2638  std::string homedir(TRAIL_CWD());
2639  std::string format("HDF");
2640  std::string newpath;
2641  int rank = 0;
2642  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
2643  std::vector<int> bcflags(3);
2644  bcflags[0] = 0;
2645  bcflags[1] = 1;
2646  bcflags[2] = 2;
2647  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocface, "RFC");
2648  int RFC_transfer = COM_get_function_handle("RFC.least_squares_transfer");
2649  int RFC_read = COM_get_function_handle( "RFC.read_overlay");
2650  if(ouf)
2651  *ouf << "Reading source window (" << trailwin << ") from "
2652  << srcpath << "/" << srcfile << "...";
2653  TRAIL_CD(srcpath,ouf);
2654  TRAIL_File2Window(srcfile,trailwin,bcflags,comm,false,true,false);
2655  newpath.assign(homedir+"/"+trgpath);
2656  MPI_Barrier(comm);
2657  if(ouf)
2658  *ouf << "done" << std::endl
2659  << "Reading target window from " << newpath
2660  << "/" << trgfile << "...";
2661  TRAIL_CD(newpath,ouf);
2662  TRAIL_File2Window(trgfile,trg,bcflags,comm,false,true,false);
2663  MPI_Barrier(comm);
2664  if(ouf)
2665  *ouf << "done" << std::endl
2666  << "Reading destination window " << r_trg << " from "
2667  << trgfile << "...";
2668  TRAIL_File2Window(trgfile,r_trg,bcflags,comm,false,true,true);
2669  MPI_Barrier(comm);
2670  newpath.assign(homedir+"/"+crpath);
2671  if(ouf)
2672  *ouf << "done" << std::endl
2673  << "Reading mesh overlay from " << newpath << "..." << std::endl;
2674  TRAIL_CD(newpath,ouf);
2675  MPI_Barrier(comm);
2676  int srcmesh = COM_get_attribute_handle( (trailwin+".mesh").c_str());
2677  int trgmesh = COM_get_attribute_handle( (trg+".mesh").c_str());
2678  std::vector<int> pane_id;
2679  if(ouf)
2680  *ouf << "TRAIL_AutoSurfer: Reading mesh overlay for all surfaces."
2681  << "TRAIL_AutoSurfer: CR DIR: " << TRAIL_CWD() << std::endl;
2682  COM_call_function( RFC_read, &srcmesh, &trgmesh, &comm,
2683  trailwin.c_str(),trg.c_str(),format.c_str());
2684  if(ouf)
2685  *ouf << "TRAIL_AutoSurfer: Beginning transfer for all surfaces."
2686  << std::endl;
2687  MPI_Barrier(comm);
2688  if(ouf)
2689  *ouf << "done" << std::endl
2690  << "Transferring data ..." << std::endl;
2691  int num_attributes;
2692  std::string names;
2693  COM_get_attributes( trailwin.c_str(),&num_attributes,names);
2694  char loc;
2695  COM_Type comtype;
2696  int ncomp;
2697  std::string unit;
2698  std::istringstream Istr(names);
2699  for(int i = 0;i < num_attributes;i++){
2700  std::string aname;
2701  Istr >> aname;
2702  COM_get_attribute(trailwin+"."+aname,&loc,&comtype,&ncomp,&unit);
2703  if((loc == 'e' || loc == 'n') && comtype == COM_DOUBLE){
2704  if(!rank)
2705  if(ouf)
2706  std::cout << "Roctrail> Transferring attribute: " << aname << " on "
2707  << (loc == 'e' ? "elements" : "nodes") << "."
2708  << std::endl;
2709  COM_resize_array((trailwin+"."+aname).c_str());
2710  COM_new_attribute((trg+"."+aname).c_str(),(char)loc,
2711  COM_DOUBLE,(int)ncomp,unit.c_str());
2712  COM_new_attribute((r_trg+"."+aname).c_str(),(char)loc,
2713  COM_DOUBLE,(int)ncomp,unit.c_str());
2714  COM_resize_array((r_trg+"."+aname).c_str());
2715  COM_resize_array((trg+"."+aname).c_str());
2716  int src_ahdl = COM_get_attribute_handle((trailwin+"."+aname).c_str());
2717  int trg_ahdl = COM_get_attribute_handle((trg+"."+aname).c_str());
2718  COM_call_function( RFC_transfer, &src_ahdl, &trg_ahdl);
2719  int *srcpane_ids;
2720  int npanes;
2721  COM_get_panes( trg.c_str(), &npanes, &srcpane_ids);
2722  pane_id.resize(npanes);
2723  for(int i = 0;i < npanes;i++)
2724  pane_id[i] = srcpane_ids[i];
2725  // These are no longer necessary as we've duped the info into
2726  // a locally allocated array
2727  COM_free_buffer( &srcpane_ids);
2728  for(int p = 0;p < npanes;p++){
2729  void *src_ptr = NULL;
2730  int src_std = 0;
2731  int src_cap = 0;
2732  void *trg_ptr = NULL;
2733  int trg_std = 0;
2734  int trg_cap = 0;
2735  COM_get_array((trg+"."+aname).c_str(),pane_id[p],
2736  &src_ptr,&src_std,&src_cap);
2737  COM_get_array((r_trg+"."+aname).c_str(),pane_id[p],
2738  &trg_ptr,&trg_std,&trg_cap);
2739  if(src_ptr && trg_ptr && (trg_std*trg_cap >= src_std*src_cap)){
2740  if(ouf)
2741  *ouf << "TRAIL_AutoSurfer: Transferred " << aname << "(" << src_std
2742  << "," << src_cap << ") to " << aname << "(" << trg_std
2743  << "," << trg_cap << ")" << std::endl;
2744  memcpy(trg_ptr,src_ptr,sizeof(double)*src_std*src_cap);
2745  }
2746  else
2747  if(ouf)
2748  *ouf << "TRAIL_AutoSurfer: WARNING: non matching sizes for "
2749  << aname << " on pane " << pane_id[p] << "."
2750  << std::endl
2751  << "TRAIL_AutoSurfer: src(" << src_std << "," << src_cap
2752  << ") trg(" << trg_std << "," << trg_cap << ")"
2753  << std::endl;
2754  }
2755  }
2756  }
2757  MPI_Barrier(comm);
2758  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocface,"RFC");
2759  newpath.assign(homedir+"/"+destpath);
2760  TRAIL_CD(newpath,ouf);
2761  MPI_Barrier(comm);
2762  if(ouf){
2763  *ouf << "Transfer is complete, massaging and writing results to "
2764  << newpath << "/" << r_trg << std::endl;
2765  if(!rank)
2766  std::cout << "Roctrail> Transfer complete, writing new surface"
2767  << std::endl;
2768  }
2769  COM_delete_window(trailwin);
2770  COM_delete_window(trg);
2771 
2772  // FIXME - Temporarily done here as a quick fix
2773  TRAIL_FixRocstarData(r_trg,ouf);
2774  TRAIL_ExtractSurf0(r_trg,"surf0",ouf);
2775  TRAIL_WriteWindow("surf0",".","surf0",".",t,id,comm);
2776  // FIXME
2777 
2778  TRAIL_WriteWindow(r_trg,".",r_trg,".",t,id,comm);
2779  TRAIL_CD(homedir,ouf);
2780  MPI_Barrier(comm);
2781  COM_delete_window(r_trg);
2782  if(ouf){
2783  *ouf << "Results written." << std::endl
2784  << "Surface data transfer is complete." << std::endl;
2785  if(!rank)
2786  std::cout << "Roctrail> New surface written" << std::endl;
2787  }
2788  MPI_Barrier(comm);
2789  return(true);
2790 }
2791 
2792 void
2793 TRAIL_WriteRocinControl(std::vector<int> &pane_id,const std::string &pre,
2794  int rank)
2795 {
2796  std::ofstream Ouf;
2797  std::string controlfilename(pre + "_in.txt");
2798  Ouf.open(controlfilename.c_str());
2799  Ouf << "@Proc: " << rank << std::endl
2800  << "@Files: " << pre << ".hdf" << std::endl;
2801  Ouf.clear();
2802  Ouf << "@Panes: ";
2803  std::vector<int>::iterator pii = pane_id.begin();
2804  while(pii != pane_id.end())
2805  Ouf << *pii++ << " ";
2806  Ouf << std::endl;
2807  Ouf.close();
2808 }
2809 
2810 bool
2811 TRAIL_WriteWindow(const std::string &wname,const std::string &winpath,
2812  const std::string &cntl_name,const std::string &cntl_path,
2813  double t,unsigned int id,MPI_Comm comm,std::ostream *ouf)
2814 {
2815 
2816  std::string timestring(TRAIL_TimeString(t));
2817  std::string homedir(TRAIL_CWD());
2818  int rank = 0;
2819  int nproc = 1;
2820  if(comm != MPI_COMM_NULL){
2821  MPI_Comm_rank(comm,&rank);
2822  MPI_Comm_size(comm,&nproc);
2823  }
2824  int *srcpane_ids;
2825  int npanes;
2826  std::vector<int> pane_id;
2827  if(ouf)
2828  *ouf << "TRAIL_WriteWindow: Entry" << std::endl;
2829  COM_get_panes( wname.c_str(), &npanes, &srcpane_ids);
2830  pane_id.resize(npanes);
2831  for(int i = 0;i < npanes;i++)
2832  pane_id[i] = srcpane_ids[i];
2833  COM_free_buffer( &srcpane_ids);
2834 
2835 
2836  COM_LOAD_MODULE_STATIC_DYNAMIC(Rocout,"Rocout");
2837  int OUT_set_option = COM_get_function_handle( "Rocout.set_option");
2838  std::string rankstr("0");
2839  COM_call_function( OUT_set_option, "rankwidth", rankstr.c_str());
2840  int whand = COM_get_function_handle("Rocout.write_attribute");
2841 
2842 
2843  int all = COM_get_attribute_handle((wname+".all"));
2844  std::ostringstream Ostr;
2845  Ostr << wname << "_" << timestring << "_" << std::setw(5)
2846  << std::setfill('0') << id;
2847  TRAIL_CD(winpath,ouf);
2848  COM_call_function(whand,Ostr.str().c_str(),&all,
2849  wname.c_str(),timestring.c_str());
2850  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocout,"Rocout");
2851 
2852  TRAIL_WriteRocinControl(pane_id,Ostr.str(),rank);
2853  if(ouf)
2854  *ouf << "TRAIL_WriteWindow: Wrote window " << wname << " id("
2855  << id << ") to " << Ostr.str() << ".hdf" << std::endl
2856  << "TRAIL_WriteWindow: Merging Rocin control files."
2857  << std::endl;
2858  if(comm != MPI_COMM_NULL)
2859  MPI_Barrier(comm);
2860  // Merge Rocin control files
2861  if(!rank)
2862  TRAIL_MergeRocinFiles(wname,wname,".",t,nproc,ouf);
2863  TRAIL_CD(homedir,ouf);
2864  if(comm != MPI_COMM_NULL)
2865  MPI_Barrier(comm);
2866  if(ouf)
2867  *ouf << "TRAIL_WriteWindow: Exit" << std::endl;
2868  return(true);
2869 }
2870 
2871 // Specify Rocout.remesh, time, and solver name
2872 double
2873 TRAIL_FindSourceTime(const std::string &dirpre,
2874  double t,
2875  const std::string &relpath)
2876 {
2877  double targ_time = -1;
2878  Directory sd(relpath);
2879  Directory::iterator dfi = sd.begin();
2880  while(dfi != sd.end()){
2881  std::string cdir(*dfi++);
2882  std::string::size_type x = cdir.find(dirpre);
2883  if(x != std::string::npos){
2884  double tt = TRAIL_TimeString(cdir.substr(dirpre.size(),9));
2885  if(tt < t && tt > targ_time)
2886  targ_time = tt;
2887 
2888  }
2889  }
2890  return(targ_time);
2891 }
2892 
2893 void
2894 TRAIL_FixRocstarData(const std::string &wname,
2895  std::ostream *ouf)
2896 {
2897  int *srcpane_ids;
2898  int npanes;
2899  std::vector<int> pane_id;
2900  COM_get_panes( wname.c_str(), &npanes, &srcpane_ids);
2901  pane_id.resize(npanes);
2902  for(int i = 0;i < npanes;i++)
2903  pane_id[i] = srcpane_ids[i];
2904  COM_free_buffer( &srcpane_ids);
2905  for(int i = 0;i < npanes;i++){
2906  double *nc_t0 = NULL;
2907  double *nc = NULL;
2908  double *mdot = NULL;
2909  int *bflag = NULL;
2910  int stride1 = 0;
2911  int stride2 = 0;
2912  int cap1 = 0;
2913  int cap2 = 0;
2914  COM_get_array((wname+".nc_t0").c_str(),pane_id[i],
2915  &nc_t0,&stride1,&cap1);
2916  COM_get_array((wname+".nc").c_str(),pane_id[i],
2917  &nc,&stride2,&cap2);
2918  if(nc_t0){
2919  if(ouf)
2920  *ouf << "TRAIL_FixRocstarData: Fixing nc_t0 for pane(" << pane_id[i]
2921  << ")" << std::endl;
2922  for(int c = 0;c < cap1;c++){
2923  unsigned int index = c*stride1;
2924  if(nc_t0[index] == 0.0 &&
2925  nc_t0[index+1] == 0.0 &&
2926  nc_t0[index+2] == 0.0){
2927  nc_t0[index] = nc[index];
2928  nc_t0[index+1] = nc[index+1];
2929  nc_t0[index+2] = nc[index+2];
2930  }
2931  }
2932  }
2933 
2934  COM_get_array((wname+".bflag").c_str(),pane_id[i],
2935  &bflag,&stride1,&cap1);
2936  COM_get_array((wname+".mdot").c_str(),pane_id[i],
2937  &mdot,&stride2,&cap2);
2938  if(!bflag && mdot){
2939  COM_resize_array((wname+".bflag").c_str());
2940  COM_get_array((wname+".bflag").c_str(),pane_id[i],
2941  &bflag,&stride1,&cap1);
2942  }
2943  if(bflag && mdot){
2944  if(ouf)
2945  *ouf << "TRAIL_FixRocstarData: Fixing bflag for pane(" << pane_id[i]
2946  << ")" << std::endl;
2947  for(int c = 0;c < cap2;c++){
2948  if(mdot[c] > 0.0)
2949  bflag[c] = 1;
2950  }
2951  }
2952  }
2953  COM_window_init_done(wname);
2954 }
2955 
2956 void
2957 TRAIL_ExtractSurf0(const std::string &srcwin,
2958  const std::string &trgwin,
2959  std::ostream *ouf)
2960 {
2961  COM_new_window(trgwin);
2962  COM_clone_attribute( (trgwin+".mesh").c_str(),
2963  (srcwin+".mesh").c_str(),1);
2964  COM_clone_attribute( (trgwin+".bcflag").c_str(),
2965  (srcwin+".bcflag").c_str(),1);
2966  int *srcpane_ids;
2967  int npanes;
2968  std::vector<int> pane_id;
2969  COM_get_panes( trgwin.c_str(), &npanes, &srcpane_ids);
2970  pane_id.resize(npanes);
2971  for(int i = 0;i < npanes;i++)
2972  pane_id[i] = srcpane_ids[i];
2973  COM_free_buffer( &srcpane_ids);
2974  for(int i = 0;i < npanes;i++){
2975  double *nc_t0 = NULL;
2976  double *nc = NULL;
2977  int stride1 = 0;
2978  int stride2 = 0;
2979  int cap1 = 0;
2980  int cap2 = 0;
2981  COM_get_array((srcwin+".nc_t0").c_str(),pane_id[i],
2982  &nc_t0,&stride1,&cap1);
2983  COM_get_array((trgwin+".nc").c_str(),pane_id[i],
2984  &nc,&stride2,&cap2);
2985  if(nc && nc_t0)
2986  memcpy(nc,nc_t0,sizeof(double)*stride2*cap2);
2987  }
2988  COM_window_init_done(trgwin);
2989 }
2990 
2991 void
2992 TRAIL_RocmopSmooth(GEM_Partition &gp, unsigned int niter)
2993 {
2994  std::string wname("smooth_win");
2995  unsigned int nnodes = gp._nc.size()/3;
2996  if(gp._debug && gp._out)
2997  *gp._out << "TRAIL_RocmopSmooth: This partition has " << nnodes
2998  << " nodes." << std::endl;
2999  MPI_Barrier(gp._comm);
3000  std::vector<double> disp(nnodes*3,0.0);
3001  if(gp._debug && gp._out)
3002  *gp._out << "TRAIL_RocmopSmooth: Creating window." << std::endl;
3003  MPI_Barrier(gp._comm);
3004  COM_new_window(wname);
3005  gp.PopulateVolumeWindow(wname);
3006  COM_new_attribute((wname+".disp"),'n',COM_DOUBLE,3,"m");
3007  COM_set_array((wname+".disp"),gp.pane_id,&disp[0]);
3008  COM_window_init_done(wname);
3009  int meshhandle = COM_get_attribute_handle((wname+".pmesh"));
3010  int disphandle = COM_get_attribute_handle((wname+".disp"));
3011  MPI_Barrier(gp._comm);
3012  if(gp._debug && gp._out)
3013  *gp._out << "TRAIL_RocmopSmooth: Loading Rocmop." << std::endl;
3014  MPI_Barrier(gp._comm);
3015  COM_LOAD_MODULE_STATIC_DYNAMIC(Rocmop,"MOP");
3016  int smoothhandle = COM_get_function_handle("MOP.smooth");
3017  int handleOption = COM_get_function_handle("MOP.set_value");
3018  MPI_Barrier(gp._comm);
3019  if(gp._debug && gp._out)
3020  *gp._out << "TRAIL_RocmopSmooth: Setting Rocmop options." << std::endl;
3021  MPI_Barrier(gp._comm);
3022  int inverted = 1;
3023  COM_call_function(handleOption,2,"inverted",&inverted) ;
3024  MPI_Barrier(gp._comm);
3025  unsigned int count = 0;
3026  while(count++ < niter){
3027  if(gp._debug && gp._out)
3028  *gp._out << "TRAIL_RocmopSmooth: Calling Rocmop's smoothing function.."
3029  << std::endl;
3030  MPI_Barrier(gp._comm);
3031  COM_call_function(smoothhandle,2,&meshhandle,&disphandle);
3032  MPI_Barrier(gp._comm);
3033  if(gp._debug && gp._out)
3034  *gp._out << "TRAIL_RocmopSmooth: Smoothing done. Applying displacements."
3035  << std::endl;
3036  unsigned int n = 0;
3037  while(n <3*nnodes){
3038  gp._nc[n] += disp[n];
3039  disp[n] = 0.0;
3040  n++;
3041  }
3042  MPI_Barrier(gp._comm);
3043  if(gp._debug && gp._out)
3044  *gp._out << "TRAIL_RocmopSmooth: Applying displacements done."
3045  << std::endl;
3046  MPI_Barrier(gp._comm);
3047  }
3048  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocmop,"MOP");
3049  COM_delete_window(wname);
3050  MPI_Barrier(gp._comm);
3051  if(gp._debug && gp._out)
3052  *gp._out << "TRAIL_RocmopSmooth: All done." << std::endl;
3053  MPI_Barrier(gp._comm);
3054 }
3055 
3056 void
3057 TRAIL_RocpropSmoothSurf(double *nc,unsigned int nnodes,
3058  unsigned int *ec,unsigned int nel,
3059  unsigned int *cnstr_type,
3060  unsigned int niter)
3061 {
3062  std::string wname("smooth_win");
3063  COM_new_window(wname,MPI_COMM_SELF);
3064  COM_new_attribute((wname+".speed"),'e',COM_DOUBLE,1,"m/s");
3065  COM_new_attribute((wname+".offsets"),'n',COM_DOUBLE,3,"m");
3066  COM_new_attribute((wname+".cnstr_type"),'e',COM_INT,1,"");
3067 
3068 
3069  for(unsigned int i = 0; i < nel*3;i++) ec[i]++;
3070  std::vector<double> offsets;
3071  offsets.resize(nnodes*3);
3072  std::vector<double> speed;
3073  speed.resize(nel,0.0);
3074 
3075  COM_set_size((wname+".nc"),102,nnodes);
3076  COM_set_array((wname+".nc"),102,nc,3);
3077  COM_set_size((wname+".:t3:real"),102,nel);
3078  COM_set_array((wname+".:t3:real"),102,ec,3);
3079  // COM_set_size((wname+".speed"),102,nel);
3080  COM_set_array((wname+".speed"),102,&speed[0],3);
3081  // COM_set_size((wname+".offsets"),102,nnodes,3);
3082  COM_set_array((wname+".offsets"),102,&offsets[0],3);
3083  // COM_set_size((wname+".cnstr_type"),102,nel);
3084  COM_set_array((wname+".cnstr_type"),102,&cnstr_type[0],1);
3085 
3086  int cnstrtype = COM_get_attribute_handle((wname+".cnstr_type"));
3087  int offset_hdl = COM_get_attribute_handle((wname+".offsets"));
3088  int speed_hdl = COM_get_attribute_handle((wname+".speed"));
3089 
3090  COM_LOAD_MODULE_STATIC_DYNAMIC( Rocprop, "PROP");
3091  int setopt = COM_get_function_handle( "PROP.set_option");
3092  COM_call_function( setopt, "method", "fo");
3093  std::ostringstream Ostr;
3094  Ostr << niter << std::ends;
3095  COM_call_function( setopt, "rediter", Ostr.str().c_str());
3096  COM_call_function( setopt, "fangle", "35");
3097  int init = COM_get_function_handle("PROP.initialize");
3098  int prop = COM_get_function_handle("PROP.propagate");
3099  int setcnstr = COM_get_function_handle("PROP.set_constraints");
3100 
3101  int pmesh_hdl = COM_get_attribute_handle((wname+".mesh"));
3102  COM_call_function(init, &pmesh_hdl);
3103 
3104  COM_call_function(setcnstr,&cnstrtype);
3105 
3106  double dt = 0.0;
3107 
3108  COM_call_function(prop,&pmesh_hdl,&speed_hdl,&dt,&offset_hdl);
3109 
3110  for(unsigned int i = 0; i < nnodes*3; i++)
3111  nc[i] += offsets[i];
3112 
3113  COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocprop,"PROP");
3114  COM_delete_window(wname);
3115  for(unsigned int i = 0; i < nel*3;i++) ec[i]--;
3116 }
3117 
3118 
Mesh::Connectivity::Sync
void Sync()
Definition: Mesh.C:246
meshes
std::vector< Mesh::UnstructuredMesh > meshes
The list of grids to write.
Definition: hdf2pltV2.C:61
COM_Type
int COM_Type
Indices for derived data types.
Definition: roccom_basic.h:122
COM_get_communicator
void COM_get_communicator(const char *wname, MPI_Comm *comm)
Definition: roccom_c++.h:346
Mesh::GenericCell_2
Definition: Mesh.H:556
TRAIL_RocmopSmooth
void TRAIL_RocmopSmooth(GEM_Partition &gp, unsigned int niter)
Definition: TRAIL.C:2992
MPI_COMM_WORLD
here we put it at the!beginning of the common block The point to point and collective!routines know about but MPI_TYPE_STRUCT as yet does not!MPI_STATUS_IGNORE and MPI_STATUSES_IGNORE are similar objects!Until the underlying MPI library implements the C version of these are declared as arrays of MPI_STATUS_SIZE!The types and are OPTIONAL!Their values are zero if they are not available Note that!using these reduces the portability of MPI_IO INTEGER MPI_BOTTOM INTEGER MPI_DOUBLE_PRECISION INTEGER MPI_LOGICAL INTEGER MPI_2REAL INTEGER MPI_2DOUBLE_COMPLEX INTEGER MPI_LB INTEGER MPI_WTIME_IS_GLOBAL INTEGER MPI_COMM_WORLD
Definition: SourceIMP/utilities/RocfracPrep/mpif.h:107
TRAIL_GetWindowSolnMetaData
void TRAIL_GetWindowSolnMetaData(const std::string &wname, std::vector< SolnMetaData > &smdv, int verblevel=0)
Definition: TRAIL.C:2260
Mesh::Connectivity::AddElement
void AddElement(const std::vector< IndexType > &elem)
Definition: Mesh.C:235
Mesh::Connectivity
General connectivity object.
Definition: Mesh.H:334
MPI_COMM_SELF
here we put it at the!beginning of the common block The point to point and collective!routines know about but MPI_TYPE_STRUCT as yet does not!MPI_STATUS_IGNORE and MPI_STATUSES_IGNORE are similar objects!Until the underlying MPI library implements the C version of these are declared as arrays of MPI_STATUS_SIZE!The types and are OPTIONAL!Their values are zero if they are not available Note that!using these reduces the portability of MPI_IO INTEGER MPI_BOTTOM INTEGER MPI_DOUBLE_PRECISION INTEGER MPI_LOGICAL INTEGER MPI_2REAL INTEGER MPI_2DOUBLE_COMPLEX INTEGER MPI_LB INTEGER MPI_WTIME_IS_GLOBAL INTEGER MPI_COMM_SELF
Definition: SourceIMP/utilities/RocfracPrep/mpif.h:107
TRAIL_HDF2Window
void TRAIL_HDF2Window(const std::string &fname, const std::string &wname, int verb=0)
Definition: TRAIL.C:2229
TRAIL_FE2FD_Transfer
int TRAIL_FE2FD_Transfer(const std::string &fewin, const std::string &fdwin, const std::string &attlist, MPI_Comm communicator, std::ostream *ouf=NULL)
Definition: TRAIL.C:974
COM_delete_window
void COM_delete_window(const char *wname)
Definition: roccom_c++.h:94
COM_INTEGER
Definition: roccom_basic.h:133
Mesh::BSExtent::CreateUnstructuredMesh
void CreateUnstructuredMesh(Connectivity &conn)
Definition: BSMesh.H:199
COM_get_attribute
void COM_get_attribute(const std::string wa_str, char *loc, int *type, int *ncomp, std::string *unit)
Definition: roccom_c++.h:269
TRAIL_FD2FE_Transfer
int TRAIL_FD2FE_Transfer(const std::string &, const std::string &, const std::string &, std::ostream *ouf=NULL)
Definition: TRAIL.C:1183
Mesh::BSExtent::FindSharedNodes
void FindSharedNodes(std::vector< BSExtent< T > > &extent_pool, std::vector< BSExtent< T > > &shared_extents, std::vector< T > &neighbors)
Definition: BSMesh.H:174
Mesh::BSExtent::Overlap
void Overlap(const BSExtent< T > &inextent, BSExtent< T > &outextent)
Definition: BSMesh.H:152
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
TRAIL_TransferSurfDataFILE
bool TRAIL_TransferSurfDataFILE(const std::string &src, const std::string &trg, const std::string &dest, const std::string &srcpath, const std::string &trgpath, const std::string &destpath, const std::string &crpath, double t, unsigned int id, MPI_Comm comm, std::ostream *=NULL)
Definition: TRAIL.C:2620
TRAIL_CreateRobustFC_old
void TRAIL_CreateRobustFC_old(const std::string &wname, const std::string &path)
Definition: TRAIL.C:2527
Mesh::BSExtent::Sync
void Sync()
Definition: BSMesh.H:88
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
COM_set_default_communicator
void COM_set_default_communicator(MPI_Comm comm)
Definition: roccom_c++.h:67
roccom.h
This file contains the prototypes for Roccom API.
TRAIL_CreateRobustFC
void TRAIL_CreateRobustFC(const std::string &wname, const std::string &path)
Definition: TRAIL.C:2515
TRAIL_AutoSurfer
void TRAIL_AutoSurfer(const std::string &src, const std::string &trg, const std::string &srcpath=".", const std::string &trgpath=".", const std::string &destpath=".", double t=0, MPI_Comm comm=MPI_COMM_NULL, std::ostream *=NULL)
Definition: TRAIL.C:2544
COM_get_array
void COM_get_array(const char *wa_str, int pane_id, void **addr, int *strd, int *cap)
Get the address for an attribute on a specific pane.
Rocmop
A Roccom mesh optimization module.
Definition: Rocmop.h:41
TRAIL_CD
int TRAIL_CD(const std::string &path, std::ostream *=NULL)
Definition: TRAIL_UnixUtils.C:181
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
COM_get_attribute_handle
int COM_get_attribute_handle(const char *waname)
Definition: roccom_c++.h:412
COM_delete_attribute
void COM_delete_attribute(const char *wa_str)
Delete an existing attribute.
Definition: roccom_c++.h:128
Mesh::GenericCell_2::ReOrient
void ReOrient(std::vector< IndexType > &ec)
Definition: Mesh.C:856
TRAIL_SetAttribute
void TRAIL_SetAttribute(const std::string &aname, int pane_id, T &value)
Definition: TRAIL.C:386
TRAIL_ExtractPanes
int TRAIL_ExtractPanes(const std::string &window_name, const std::string &attribute_name, int attribute_value, std::vector< int > &pane_id)
Definition: TRAIL.C:129
TRAIL_RocpropSmoothSurf
void TRAIL_RocpropSmoothSurf(double *nc, unsigned int nnodes, unsigned int *ec, unsigned int nel, unsigned int *cnstr_type, unsigned int niter)
Definition: TRAIL.C:3057
smdv
std::vector< SolnMetaData > smdv
The solution Metadata.
Definition: hdf2pltV2.C:63
TRAIL_FD2FE_WinCreate
int TRAIL_FD2FE_WinCreate(const std::string &wname, const std::string &outwname, std::ostream *ouf=NULL)
Adds ghost zones for block structured meshes to close gaps in the interface surface mesh...
Definition: TRAIL.C:565
TRAIL_Window2UnstructuredMesh
void TRAIL_Window2UnstructuredMesh(const std::string &wname, std::vector< Mesh::UnstructuredMesh > &meshes, std::vector< SolnMetaData > &smdv, std::vector< std::vector< std::vector< double > > > &soln_data, int verb=0, bool no_ghost=false)
Definition: TRAIL.C:2010
TRAIL_TimeString
std::string TRAIL_TimeString(double t)
Definition: TRAIL.C:85
copy
*********************************************************************Illinois Open Source License ****University of Illinois NCSA **Open Source License University of Illinois All rights reserved ****Developed free of to any person **obtaining a copy of this software and associated documentation to deal with the Software without including without limitation the rights to ** copy
Definition: roccomf90.h:20
TRAIL_FixRocstarData
void TRAIL_FixRocstarData(const std::string &wname, std::ostream *ouf=NULL)
Definition: TRAIL.C:2894
COM_UNLOAD_MODULE_STATIC_DYNAMIC
#define COM_UNLOAD_MODULE_STATIC_DYNAMIC(moduleName, windowString)
Definition: roccom_basic.h:113
Mesh::BSExtent::Flatten
void Flatten(std::vector< T > &output)
Definition: BSMesh.H:111
TRAIL_UniqueAcrossProcs
int TRAIL_UniqueAcrossProcs(std::vector< int > &input_data, std::vector< int > &output_data, MPI_Comm communicator)
Definition: TRAIL.C:191
Mesh::NodalCoordinates::Size
IndexType Size() const
Definition: Mesh.C:49
Mesh::BSExtent::GetFlatIndices
void GetFlatIndices(const BSExtent< T > &extent, std::vector< T > &indices)
Definition: BSMesh.H:139
COM_get_attributes
void COM_get_attributes(const char *wname, int *na, std::string &names)
Definition: roccom_c++.h:360
GEM_Partition
Definition: GEM.H:282
GEM_Partition::_id
unsigned int _id
Definition: GEM.H:284
GEM_Partition::_out
std::ostream * _out
Definition: GEM.H:308
TRAIL_CWD
std::string TRAIL_CWD(void)
Definition: TRAIL_UnixUtils.C:174
TRAIL_GetRocstarDumpStrings
void TRAIL_GetRocstarDumpStrings(const std::string &filename, std::string &wname, std::string &timestring, std::string &rankstring)
Definition: TRAIL.C:112
GEM_Partition::pane_id
int pane_id
Definition: GEM.H:315
Rocin
Definition: Rocin.h:64
Rocout
Definition: Rocout.h:81
TRAIL_Search_Block_Structured_Pool
int TRAIL_Search_Block_Structured_Pool(std::vector< std::vector< int > > &search_extent, std::vector< std::vector< std::vector< int > > > &extent_pool, std::vector< std::vector< std::vector< int > > > &neighbor_extent, std::vector< int > &neighbors)
Given an array of ranges: search_extent[nd][2] and one array of array of ranges: partition_extent[num...
Definition: TRAIL.C:229
i
blockLoc i
Definition: read.cpp:79
TRAIL_Att2Vec
int TRAIL_Att2Vec(const std::string &att, int pane_id, std::vector< T > &dest)
Definition: TRAIL.C:952
x
void int int REAL * x
Definition: read.cpp:74
Mesh::Connectivity::Nelem
IndexType Nelem() const
Definition: Mesh.H:364
COM_window_init_done
void COM_window_init_done(const char *w_str, int pane_changed=true)
Definition: roccom_c++.h:102
TRAIL_UnstructuredMesh2Pane
int TRAIL_UnstructuredMesh2Pane(const std::string &wname, int pane_id, Mesh::UnstructuredMesh &mesh, SolnMetaData &smdv, std::vector< std::vector< double > > &soln_data, int verblevel)
Creates a window from a Mesh object. (copies data)
Definition: TRAIL.C:473
n
const NT & n
Definition: rational_rotation.h:72
COM_get_size
void COM_get_size(const char *wa_str, int pane_id, int *size, int *ng=0)
Get the sizes of an attribute.
Definition: roccom_c++.h:274
COM_get_sizeof
int COM_get_sizeof(const COM_Type type, int c)
Definition: roccom_c++.h:560
COM_clone_attribute
void COM_clone_attribute(const char *wname, const char *attr, int wg=1, const char *ptnname=0, int val=0)
Clone the subset of panes of another window of which the given pane attribute has value val...
Definition: roccom_c++.h:234
COM_new_window
void COM_new_window(const char *wname, MPI_Comm c=MPI_COMM_NULL)
Definition: roccom_c++.h:86
COM_call_function
void COM_call_function(const int wf, int argc,...)
Definition: roccom_c.C:48
TRAIL_SurfaceMesh2Window
int TRAIL_SurfaceMesh2Window(const std::string &wname, int pane_id, Mesh::NodalCoordinates &, Mesh::Connectivity &)
Creates a window from a Mesh object. (copies data)
Definition: TRAIL.C:414
GEM_Partition::debug
bool debug(bool s=true)
Definition: GEM.C:1266
TRAIL_FindSourceTime
double TRAIL_FindSourceTime(const std::string &dirpre, double t, const std::string &relpath)
Definition: TRAIL.C:2873
TRAIL_File2Window
void TRAIL_File2Window(const std::string &fname, const std::string &wname, std::vector< int > &bcflags, MPI_Comm comm=MPI_COMM_NULL, bool apply_disp=false, bool all=false, bool with_ghost=false)
Definition: TRAIL.C:2391
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
COM_INT
Definition: roccom_basic.h:130
Mesh::Connectivity::SyncSizes
void SyncSizes()
Definition: Mesh.C:281
TRAIL_CreateDirectory
int TRAIL_CreateDirectory(const std::string &fname)
Definition: TRAIL_UnixUtils.C:81
Mesh::Connectivity::Flatten
void Flatten(std::vector< T > &outcon) const
Definition: Mesh.H:394
min
Vector_n min(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:346
TRAIL_Get_Block_Structured_Neighbors
int TRAIL_Get_Block_Structured_Neighbors(std::vector< std::vector< int > > &local_extent, std::vector< std::vector< int > > &global_extent, std::vector< std::vector< std::vector< int > > > &extent_pool, std::vector< std::vector< int > > &ghost_extent, std::vector< std::vector< std::vector< int > > > &neighbor_extent, std::vector< int > &neighbors)
Given two arrays of ranges: global_extent[nd][2] local_extent[nd][2] an extent pool: extent_pool[npoo...
Definition: TRAIL.C:277
TRAIL_Add_Attributes
int TRAIL_Add_Attributes(const std::string &srcwin, const std::string &destwin, const std::vector< std::string > &atts)
Definition: TRAIL.C:1156
COM_delete_pane
void COM_delete_pane(const char *str, int pid)
Definition: roccom_c++.h:110
j
j indices j
Definition: Indexing.h:6
id
unsigned long id(const Leda_like_handle &x)
Definition: Handle.h:107
pi
const double pi
Definition: compare_normals.C:253
Rocmap
Definition: Rocmap.h:39
GEM_Partition::_comm
MPI_Comm _comm
Definition: GEM.H:312
GEM_Partition::_debug
bool _debug
Definition: GEM.H:310
TRAIL_MergeRocinFiles
void TRAIL_MergeRocinFiles(const std::string &srcname, const std::string &trgname, const std::string &path="./", double t=0, unsigned int np=1, std::ostream *=NULL)
Definition: TRAIL.C:2476
TRAIL_Debug
void TRAIL_Debug(GEM_Partition &gp)
Definition: TRAIL.C:60
Mesh::UnstructuredMesh::con
Connectivity con
Definition: Mesh.H:450
GEM_Partition::PopulateVolumeWindow
bool PopulateVolumeWindow(const std::string &wname)
Definition: GEM.C:747
ni
void int int REAL REAL REAL *z blockDim dim * ni
Definition: read.cpp:77
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
COM_get_panes
void COM_get_panes(const char *wname, std::vector< int > &pane_ids, int rank=-2)
Definition: roccom_c++.h:350
COM_free_buffer
void COM_free_buffer(int **buf)
Definition: roccom_c++.h:397
rank
static int rank
Definition: advectest.C:66
TRAIL_WriteRocinControl
void TRAIL_WriteRocinControl(std::vector< int > &pane_id, const std::string &pre, int rank)
Definition: TRAIL.C:2793
Mesh::IndexType
IRAD::Primitive::IndexType IndexType
Definition: Mesh.H:57
TRAIL_GetPanelAttribute
int TRAIL_GetPanelAttribute(const std::string &window_name, const std::string &attribute_name, const std::string &qual_name, int qualval, std::vector< int > &attvec)
Definition: TRAIL.C:149
GEM_Partition::_nc
std::vector< double > _nc
Definition: GEM.H:293
COM_get_default_communicator
MPI_Comm COM_get_default_communicator()
Definition: roccom_c++.h:69
TRAIL_ExtractSurf0
void TRAIL_ExtractSurf0(const std::string &srcwin, const std::string &trgwin, std::ostream *ouf=NULL)
Definition: TRAIL.C:2957
COM_DOUBLE
Definition: roccom_basic.h:131
Mesh::BSExtent
Simple Block Structured Mesh object.
Definition: BSMesh.H:16
Mesquite::pow
double pow(double value, const Exponent &exp)
Definition: includeLinks/Exponent.hpp:89
TRAIL_Copy2Attribute
void TRAIL_Copy2Attribute(const std::string &aname, const std::vector< T > &container, int pane_id, int asize=1)
Definition: TRAIL.C:373
COM_LOAD_MODULE_STATIC_DYNAMIC
#define COM_LOAD_MODULE_STATIC_DYNAMIC(moduleName, windowString)
Definition: roccom_basic.h:111
COM_resize_array
void COM_resize_array(const char *wa_str, int pane_id=0, void **addr=NULL, int strd=-1, int cap=0)
Resize an attribute on a specific pane and return the address by setting addr.
Definition: roccom_c++.h:200
COM_get_function_handle
int COM_get_function_handle(const char *wfname)
Definition: roccom_c++.h:428
Mesh::UnstructuredMesh::nc
NodalCoordinates nc
Definition: Mesh.H:449
TRAIL_FD2FE_WinCreate2
int TRAIL_FD2FE_WinCreate2(const std::string &wname, const std::string &outwname, std::ostream *ouf=NULL)
Takes as input a block structured FD grid.
Definition: TRAIL.C:1652
TRAIL_GetWindowSolnData
void TRAIL_GetWindowSolnData(const std::string &wname, std::vector< std::vector< std::vector< double > > > &soln_data, std::vector< SolnMetaData > &smdv, int verblevel=0)
Definition: TRAIL.C:2305
COM_EXTERN_MODULE
#define COM_EXTERN_MODULE(moduleName)
Definition: roccom_basic.h:116