Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Functions
flotsam.C File Reference
#include <iostream>
#include <iomanip>
#include <limits>
#include "Global.H"
#include "PMesh.H"
#include "BSMesh.H"
#include "Profiler.H"
#include "FloGrid.H"
#include "flotsam_cl.H"
Include dependency graph for flotsam.C:

Go to the source code of this file.

Functions

void test2 (void)
 
void test (void)
 
void writeVtkData (Mesh::NodalCoordinates &nc, Mesh::Connectivity &con, const std::string &filename, std::vector< double > &soln, std::vector< Mesh::IndexType > indices, unsigned int stride)
 
int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 274 of file flotsam.C.

References FloGrid::BlockCount(), FloGrid::Blocks(), FloGrid::CloseSolutionFile(), Mesh::CollideCellsWithBox(), Mesh::CollideMeshWithBox(), FloGrid::CreateSolutions(), d, dist(), CBox::empty(), Mesh::FindPointInCells(), Mesh::BSExtent< T >::GetFlatIndices(), Mesh::GetMeshBoxes(), Mesh::GetMeshCentroids(), i, Mesh::BSExtent< T >::Init(), Mesh::NodalCoordinates::init(), FlotsamComLine::Initialize(), intersection(), Mesh::GenericElement::Inverted(), CVector::mag(), max(), CBox::merge(), nvc::norm(), FloGrid::OpenSolutionFile(), CBox::P1(), CBox::P2(), rank, FloGrid::ReadAllBlocks(), FloGrid::ReadAllSolutions(), FloGrid::ReadNBlocks(), Mesh::GenericElement::ReOrient(), FloGrid::SetGhostLayers(), FloGrid::SolnFile(), FloGrid::Time(), cimg_library::cimg::time(), FloGrid::UnknownNumber(), FloGrid::WriteBlocks(), CPoint::x(), CPoint::y(), and CPoint::z().

275 {
276  // test2();
277  // return(0);
278  std::ostream* StdOut = NULL; // program output
279  std::ostream* ErrOut = NULL; // program errors
280  std::ostream* LogOut = NULL; // log for each proc (if debugging on)
281  std::ofstream LogFile;
282  IRAD::Global::GlobalObj<std::string,Mesh::IndexType,IRAD::Profiler::ProfilerObj> global;
283  IRAD::Comm::CommunicatorObject comm(&argc,&argv);
284  unsigned int rank = comm.Rank();
285  unsigned int nproc = comm.Size();
286  // std::cout << "rank/nproc: " << rank << "/" << nproc << std::endl;
287  if(rank == 0){
288  StdOut = &std::cout;
289  ErrOut = &std::cerr;
290  }
291  global.Init("flotsam",rank);
292  global.SetDebugLevel(0);
293  global.FunctionEntry("main");
294 
295  FlotsamComLine comline((const char **)argv);
296  comline.Initialize();
297  int clerr = comline.ProcessOptions();
298  std::string sverb = comline.GetOption("verb");
299  bool debug = !comline.GetOption("debug").empty();
300  bool array_checking = !comline.GetOption("checking").empty();
301 
302  IRAD::Comm::DataTypes IndexDT = (sizeof(Mesh::IndexType) == sizeof(size_t) ?
303  IRAD::Comm::DTSIZET : IRAD::Comm::DTUINT);
304  if(!comline.GetOption("help").empty()){
305  if(StdOut)
306  *StdOut << comline.LongUsage() << std::endl;
307  comm.SetExit(1);
308  }
309  if(comm.Check())
310  return(0);
311 
312  if(clerr){
313  if(ErrOut)
314  *ErrOut << comline.ErrorReport() << std::endl
315  << std::endl << comline.ShortUsage() << std::endl;
316  comm.SetExit(1);
317  }
318  if(comm.Check())
319  return(1);
320  int verblevel = 0;
321  if(!sverb.empty()){
322  verblevel = 1;
323  if(sverb != ".true."){
324  std::istringstream Istr(sverb);
325  Istr >> verblevel;
326  }
327  }
328  if(nproc > 1){
329  if(rank == 0 && verblevel > 0)
330  std::cout << "Flotsam running on " << nproc << " processors." << std::endl;
331  std::ostringstream debugfilename;
332  debugfilename << comline.ProgramName() << ".output." << rank;
333  LogFile.open(debugfilename.str().c_str());
334  global.SetDebugStream(LogFile);
335  LogOut = &LogFile;
336  ErrOut = &LogFile;
337  if(verblevel > 1 && rank > 0)
338  StdOut = &LogFile;
339  }
340  else{
341  global.SetDebugStream(std::cerr);
342  LogOut=&std::cout;
343  ErrOut=&std::cout;
344  }
345  if(verblevel > 1 || debug){
346  global.SetDebugLevel(1);
347  }
348  if(debug && rank > 0)
349  StdOut = LogOut;
350  std::vector<std::string> infiles = comline.GetArgs();
351  if(infiles.size()==0) {
352  if(ErrOut)
353  *ErrOut << "flotsam::Error: No input specified." << std::endl
354  << std::endl << comline.ShortUsage() << std::endl;
355  comm.SetExit(1);
356  }
357  if(comm.Check())
358  return(1);
359 
360  // XXXXXXXXXXXXXXXXXXXXXXXXXX
361  // Do not modify above here.
362  // XXXXXXXXXXXXXXXXXXXXXXXXXX
363 
364  std::string casename(infiles[0]);
365  std::string timestamp(infiles[1]);
366  std::string targname(infiles[2]);
367 
368  FloGrid sourcegrid(casename);
369  FloGrid targetgrid(targname);
370 
371  unsigned int total_source_cells = 0;
372  unsigned int total_source_nodes = 0;
373  unsigned int total_target_cells = 0;
374  unsigned int total_target_nodes = 0;
375  unsigned int NSourceRegions = 0;
376  unsigned int NTargetRegions = 0;
377  unsigned int totalcells = 0;
378 
379  // This function reads the entire Rocflo grid from the file
380  // casename.grda_0.00000E-00
381  // The grid blocks are kept separated as they are in the input
382  std::vector<std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> > > target_source_cells;
383  std::vector<std::vector<Mesh::IndexType> > source_cell_ids;
384  std::vector<std::vector<Mesh::IndexType> > target_cell_ids;
385 
386  std::vector<Mesh::Connectivity> source_mesh;
387  std::vector<Mesh::Connectivity> target_mesh;
388 
389  bool do_search = true;
390  if(do_search) {
391  global.FunctionEntry("ReadFloGrid");
392  sourcegrid.ReadAllBlocks();
393  global.FunctionExit("ReadFloGrid");
394 
395 
396  // The number of regions is the number of blocks in the input grid
397  NSourceRegions = sourcegrid.Blocks().size();
398 
399  std::vector<std::vector<double> > source_cell_centers(NSourceRegions);
400 
401  // The extents are just an abstraction to represent indexing into a logically
402  // rectangular dataset ijk.
403  std::vector<Mesh::BSExtent<Mesh::IndexType> > source_extents(NSourceRegions);
404 
405  // The connectivities are the unstructured-like representation of the input
406  // block structured mesh.
407  source_mesh.resize(NSourceRegions);
408 
409  // The boxes are: boxes[region_index][0] = Box around entire region
410  // boxes[region_index][1] = " " smallest element in region
411  // boxes[region_index][2] = " " largest element in region
412  GeoPrim::CBox sourcebox;
413  std::vector<std::vector<GeoPrim::CBox> > source_boxes(NSourceRegions);
414 
415 
416 
417  // This section loops through and populates the unstructured representations and
418  // block bounding boxes.
419  std::vector<FloGridBlock>::iterator source_block_iterator = sourcegrid.Blocks().begin();
420  std::vector<Mesh::BSExtent<Mesh::IndexType> >::iterator source_extent_iterator = source_extents.begin();
421  std::vector<Mesh::Connectivity>::iterator source_mesh_iterator = source_mesh.begin();
422  std::vector<std::vector<GeoPrim::CBox> >::iterator source_boxes_iterator = source_boxes.begin();
423  while( source_block_iterator != sourcegrid.Blocks().end())
424  {
425  // Convert the block coordinates from
426  // block format to sequential format, i.e.
427  // [--x--,--y--,--z--] to [xyzxyzxyz]
428  source_block_iterator->MeshLayout();
429 
430  // Block extents go from 1 to (size+1). Size indicates
431  // number of *elements* in a given direction while
432  // the number of nodes in that direction is size+1.
433  std::vector<Mesh::IndexType> source_block_extent(6,0);
434  unsigned int isize = source_block_iterator->isize();
435  unsigned int jsize = source_block_iterator->jsize();
436  unsigned int ksize = source_block_iterator->ksize();
437  source_block_extent[0] = source_block_extent[2] = source_block_extent[4] = 1;
438  source_block_extent[1] = isize+1;
439  source_block_extent[3] = jsize+1;
440  source_block_extent[5] = ksize+1;
441  source_extent_iterator->Init(source_block_extent);
442 
443  // Use the extent object's canned routine for converting
444  // to an unstructured connectivity array
445  global.FunctionEntry("CreateUnsMesh");
446  source_extent_iterator->CreateUnstructuredMesh(*source_mesh_iterator);
447  source_mesh_iterator->Sync();
448  source_mesh_iterator->SyncSizes();
449  global.FunctionExit("CreateUnsMesh");
450 
451 
452  // Create a lightweight NodalCoordinates object around
453  // the existing block coordinates in place. This will
454  // be used in the canned unstructured routines.
455  Mesh::NodalCoordinates nc;
456  unsigned int nnodes = source_block_iterator->NNodes();
457  nc.init(nnodes,&((source_block_iterator->Coords())[0]));
458  unsigned int source_region_index = source_block_iterator - sourcegrid.Blocks().begin();
459  Mesh::GetMeshCentroids(nc,*source_mesh_iterator,source_cell_centers[source_region_index]);
460  global.FunctionEntry("Orient elements");
461  Mesh::IndexType invertedcount = 0;
462  unsigned int number_of_elements = isize*jsize*ksize;
463  for(Mesh::IndexType element_being_processed = 1;
464  element_being_processed <= number_of_elements;
465  element_being_processed++)
466  {
467  Mesh::IndexType index = element_being_processed - 1;
468  Mesh::IndexType size_of_element = (*source_mesh_iterator)[index].size();
469  Mesh::GenericElement ge(size_of_element);
470  if(ge.Inverted((*source_mesh_iterator)[index],nc)){
471  invertedcount++;
472  ge.ReOrient((*source_mesh_iterator)[index]);
473  }
474  }
475  global.FunctionExit("Orient elements");
476  // if(StdOut && verblevel)
477  // *StdOut << "Number of elements reoriented: " << invertedcount
478  // << std::endl;
479 
480  global.FunctionEntry("BoxProcessing");
481 
482  // Get the characteristic boxes on the new unstructured mesh
483  source_boxes_iterator->resize(3);
484  Mesh::GetMeshBoxes(nc,*source_mesh_iterator,(*source_boxes_iterator)[0],(*source_boxes_iterator)[1],
485  (*source_boxes_iterator)[2]);
486  sourcebox.merge((*source_boxes_iterator)[0]); // for debugging - ensure sourcebox = region box
487  global.FunctionExit("BoxProcessing");
488 
489  // track the total number of nodes and cells
490  total_source_nodes += nnodes;
491  total_source_cells += isize*jsize*ksize;
492 
493  source_boxes_iterator++;
494  source_mesh_iterator++;
495  source_extent_iterator++;
496  source_block_iterator++;
497  }
498 
499  if(StdOut && verblevel > 0)
500  *StdOut << "Source grid has " << NSourceRegions << " regions, "
501  << total_source_nodes << " nodes, and " << total_source_cells
502  << " cells." << std::endl
503  << "With bounding box: " << sourcebox << std::endl;
504 
505  //
506  // Repeat everything above on the target mesh
507  //
508  // FloGrid targetgrid(targname);
509  GeoPrim::CBox targetbox;
510  if(!rank)
511  NTargetRegions = targetgrid.BlockCount();
512  comm.BroadCast(NTargetRegions,0);
513  if(debug)
514  std::cout << "ntarget regions = " << NTargetRegions << std::endl;
515  int ntargregpproc = NTargetRegions/nproc;
516  int nleftovers = NTargetRegions%nproc;
517  int first = rank * ntargregpproc + 1;
518  int ntargreg = ntargregpproc;
519  if(rank < nleftovers){
520  ntargreg++;
521  first+=rank;
522  }
523  else
524  first += nleftovers;
525  if(debug){
526  std::cout << "First/N: " << first << "/" << ntargreg << std::endl;
527  }
528  global.FunctionEntry("ReadFloGrid");
529  // targetgrid.ReadAllBlocks();
530  targetgrid.ReadNBlocks(first,ntargreg);
531  global.FunctionExit("ReadFloGrid");
532  unsigned int TargetGridSize = NTargetRegions;
533  NTargetRegions = ntargreg;
534  // NTargetRegions = targetgrid.Blocks().size();
535  std::vector<Mesh::BSExtent<Mesh::IndexType> > target_extents(NTargetRegions);
536  target_mesh.resize(NTargetRegions);
537 
538  std::vector<std::vector<GeoPrim::CBox> > target_boxes(NTargetRegions);
539 
540  std::vector<FloGridBlock>::iterator target_block_iterator = targetgrid.Blocks().begin();
541  std::vector<Mesh::BSExtent<Mesh::IndexType> >::iterator target_extent_iterator = target_extents.begin();
542  std::vector<Mesh::Connectivity>::iterator target_mesh_iterator = target_mesh.begin();
543  std::vector<std::vector<GeoPrim::CBox> >::iterator target_boxes_iterator = target_boxes.begin();
544  while( target_block_iterator != targetgrid.Blocks().end())
545  {
546  target_block_iterator->MeshLayout();
547  std::vector<Mesh::IndexType> target_block_extent(6,0);
548  unsigned int isize = target_block_iterator->isize();
549  unsigned int jsize = target_block_iterator->jsize();
550  unsigned int ksize = target_block_iterator->ksize();
551  target_block_extent[0] = target_block_extent[2] = target_block_extent[4] = 1;
552  target_block_extent[1] = isize+1;
553  target_block_extent[3] = jsize+1;
554  target_block_extent[5] = ksize+1;
555  target_extent_iterator->Init(target_block_extent);
556  global.FunctionEntry("CreateUnsMesh");
557  target_extent_iterator->CreateUnstructuredMesh(*target_mesh_iterator);
558  target_mesh_iterator->Sync();
559  target_mesh_iterator->SyncSizes();
560  global.FunctionExit("CreateUnsMesh");
561  target_boxes_iterator->resize(3);
562  Mesh::NodalCoordinates nc;
563  unsigned int nnodes = target_block_iterator->NNodes();
564  total_target_nodes += nnodes;
565  total_target_cells += (isize*jsize*ksize);
566  nc.init(nnodes,&((target_block_iterator->Coords())[0]));
567  global.FunctionEntry("Orient elements");
568  Mesh::IndexType invertedcount = 0;
569  unsigned int number_of_elements = isize*jsize*ksize;
570  for(Mesh::IndexType element_being_processed = 1;
571  element_being_processed <= number_of_elements;
572  element_being_processed++)
573  {
574  Mesh::IndexType index = element_being_processed - 1;
575  Mesh::IndexType size_of_element = (*target_mesh_iterator)[index].size();
576  Mesh::GenericElement ge(size_of_element);
577  if(ge.Inverted((*target_mesh_iterator)[index],nc)){
578  invertedcount++;
579  ge.ReOrient((*target_mesh_iterator)[index]);
580  }
581  }
582  global.FunctionExit("Orient elements");
583  // if(StdOut && verblevel)
584  // *StdOut << "Number of elements reoriented: " << invertedcount
585  // << std::endl;
586  global.FunctionEntry("BoxProcessing");
587  Mesh::GetMeshBoxes(nc,*target_mesh_iterator,(*target_boxes_iterator)[0],(*target_boxes_iterator)[1],
588  (*target_boxes_iterator)[2]);
589  targetbox.merge((*target_boxes_iterator)[0]);
590  global.FunctionExit("BoxProcessing");
591  target_boxes_iterator++;
592  target_mesh_iterator++;
593  target_extent_iterator++;
594  target_block_iterator++;
595  }
596 
597  unsigned int totalnodes = 0;
598  unsigned int nregions = 0;
599  double minx, maxx, miny, maxy, minz, maxz;
600  minx = miny = minz = std::numeric_limits<double>::max();
601  maxx = maxy = maxz = -1 * std::numeric_limits<double>::max();
602  comm.Reduce(NTargetRegions,nregions,IRAD::Comm::DTUINT,IRAD::Comm::SUMOP,0);
603  comm.Reduce(total_target_cells,totalcells,IRAD::Comm::DTUINT,IRAD::Comm::SUMOP,0);
604  comm.Reduce(total_target_nodes,totalnodes,IRAD::Comm::DTUINT,IRAD::Comm::SUMOP,0);
605  comm.Reduce(targetbox.P1().x(),minx,IRAD::Comm::DTDOUBLE,IRAD::Comm::MINOP,0);
606  comm.Reduce(targetbox.P1().y(),miny,IRAD::Comm::DTDOUBLE,IRAD::Comm::MINOP,0);
607  comm.Reduce(targetbox.P1().z(),minz,IRAD::Comm::DTDOUBLE,IRAD::Comm::MINOP,0);
608  comm.Reduce(targetbox.P2().x(),maxx,IRAD::Comm::DTDOUBLE,IRAD::Comm::MAXOP,0);
609  comm.Reduce(targetbox.P2().y(),maxy,IRAD::Comm::DTDOUBLE,IRAD::Comm::MAXOP,0);
610  comm.Reduce(targetbox.P2().z(),maxz,IRAD::Comm::DTDOUBLE,IRAD::Comm::MAXOP,0);
611  GeoPrim::CBox domainbox(GeoPrim::CPoint(minx,miny,minz),GeoPrim::CPoint(maxx,maxy,maxz));
612 
613  if(StdOut && verblevel > 0) {
614  if(rank)
615  *StdOut << "Target grid (" << rank << ") has " << NTargetRegions << " regions, "
616  << total_target_nodes << " nodes, and " << total_target_cells
617  << " cells." << std::endl
618  << "With bounding box: " << targetbox << std::endl;
619  else
620  *StdOut << "Target grid has " << nregions << " regions, "
621  << totalnodes << " nodes, and " << totalcells
622  << " cells." << std::endl
623  << "With bounding box: " << domainbox << std::endl;
624  }
625 
626 
627  // for every target region, which source regions
628  // target_connectivity[target_region_index][1:nintersections] = source_region_id
629  Mesh::Connectivity target_connectivity(NTargetRegions);
630 
631  // holds the intesections of the colliding target/source regions
632  // target_intersections[target_region_index][1:nintersecting_regions] = intersecting_box
633  std::vector< std::vector<GeoPrim::CBox> > target_intersections(NTargetRegions);
634 
635 
636  Mesh::Connectivity::iterator target_connectivity_iterator = target_connectivity.begin();
637  std::vector< std::vector<GeoPrim::CBox> >::iterator target_intersections_iterator = target_intersections.begin();
638  target_boxes_iterator = target_boxes.begin();
639  unsigned int target_id = 0;
640 
641  // This loop collides the target and source regions, thereby populating the
642  // target connectivity and target_intersections
643  global.FunctionEntry("RegionIntersection");
644  while(target_boxes_iterator != target_boxes.end()){
645  target_id++;
646  source_boxes_iterator = source_boxes.begin();
647  unsigned int source_id = 0;
648  // std::cout << "Target region " << target_id << " collides with source regions: ";
649  while(source_boxes_iterator != source_boxes.end()){
650  source_id++;
651  GeoPrim::CBox intersection = (*target_boxes_iterator)[0].intersect((*source_boxes_iterator)[0]);
652  if(!intersection.empty()){
653  // std::cout << source_id << " ";
654  target_connectivity_iterator->push_back(source_id);
655  target_intersections_iterator->push_back(intersection);
656  }
657  source_boxes_iterator++;
658  }
659  // std::cout << std::endl;
660  target_connectivity_iterator++;
661  target_boxes_iterator++;
662  target_intersections_iterator++;
663  }
664  global.FunctionExit("RegionIntersection");
665 
666 
667  global.FunctionEntry("FindIntersectionCells");
668  // Holds the nodes that participate in each intersecting region
669  // target_intersection_nodes[1:ntargetregions][1:nintersection][1:nnodes_intersecting_region]
670  std::vector<Mesh::Connectivity> target_intersection_cells(NTargetRegions);
671  std::vector<std::vector<double> > target_region_cell_centers(NTargetRegions);
672 
673  // For every target region
674  // For every target region's cell
675  // a vector of pairs<source region id,cell id>
676  // where the target cell is found
677  std::vector<std::vector<std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> > > > target_source_matches(NTargetRegions);
678  target_source_cells.resize(NTargetRegions);
679  std::vector<Mesh::Connectivity>::iterator tint_cells_iterator = target_intersection_cells.begin();
680  std::vector<std::vector<std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> > > >::iterator tsmit =
681  target_source_matches.begin();
682  target_connectivity_iterator = target_connectivity.begin();
683  target_boxes_iterator = target_boxes.begin();
684  target_block_iterator = targetgrid.Blocks().begin();
685  target_mesh_iterator = target_mesh.begin();
686  target_intersections_iterator = target_intersections.begin();
687  while(tint_cells_iterator != target_intersection_cells.end())
688  {
689  // std::ofstream Ouf;
690  unsigned int target_region_index = tint_cells_iterator - target_intersection_cells.begin();
691  unsigned int nnodes = target_block_iterator->NNodes();
692  unsigned int ncells = target_block_iterator->NumCells();
693  // std::ostringstream Ostr;
694  // Ostr << "debug" << target_region_index+1 ;
695  // Ouf.open(Ostr.str().c_str());
696 
697  tsmit->resize(ncells);
698  target_source_cells[target_region_index].resize(ncells);
699  // std::vector<unsigned int> nfound(ncells,0);
700  std::vector<bool> cell_found(ncells,false);
701  Mesh::NodalCoordinates nc(nnodes,&(target_block_iterator->Coords()[0]));
702  GeoPrim::CBox intersection_box;
703  global.FunctionEntry("Centroids");
704  Mesh::GetMeshCentroids(nc,*target_mesh_iterator,target_region_cell_centers[target_region_index]);
705  global.FunctionExit("Centroids");
706  tint_cells_iterator->Resize(target_intersections_iterator->size());
707  Mesh::Connectivity::iterator intersection_cells_iterator = tint_cells_iterator->begin();
708  std::vector<GeoPrim::CBox>::iterator intersections_iterator = target_intersections_iterator->begin();
709  std::vector<Mesh::IndexType>::iterator source_region_ids = target_connectivity_iterator->begin();
710  while(intersections_iterator != target_intersections_iterator->end())
711  {
712  // Initialize source mesh data structures
713  Mesh::IndexType source_region_index = *source_region_ids++ - 1;
714  Mesh::IndexType nsource_nodes = sourcegrid.Blocks()[source_region_index].NNodes();
715  Mesh::NodalCoordinates source_coords(nsource_nodes,&(sourcegrid.Blocks()[source_region_index].Coords()[0]));
716  Mesh::Connectivity &sourcemesh = source_mesh[source_region_index];
717  double source_edgelen = (source_boxes[source_region_index][1].P2() - source_boxes[source_region_index][1].P1()).norm();
718  // Obtain the source cells which live in this intersection
719  std::vector<Mesh::IndexType> source_cells;
720  global.FunctionEntry("CollideMeshBox");
721  Mesh::CollideMeshWithBox(source_coords,sourcemesh,*intersections_iterator,source_cells);
722  global.FunctionExit("CollideMeshBox");
723  std::vector<Mesh::IndexType> target_cells;
724  global.FunctionEntry("CollideMeshBox");
725  Mesh::CollideMeshWithBox(nc,*target_mesh_iterator,*intersections_iterator,target_cells);
726  global.FunctionExit("CollideMeshBox");
727 
728  if(source_cells.empty() || target_cells.empty()){
729  // std::cout << "found " << source_cells.size() << " source candidates." << std::endl
730  // << "found " << target_cells.size() << " target candidates." << std::endl
731  // << "source id = " << source_region_index+1 << ", box: " << *intersections_iterator
732  // << std::endl << "source box = " << source_boxes[source_region_index][0] << std::endl;
733  }
734  else {
735  if(!rank)
736  std::cout << ".";
737  std::vector<Mesh::IndexType>::iterator ttccii = target_cells.begin();
738  // for(int i = 0;i < ncells;i++){
739  while(ttccii != target_cells.end()){
740  int i = *ttccii++ - 1;
741  intersection_box.merge(*intersections_iterator);
742  GeoPrim::CPoint p(&(target_region_cell_centers[target_region_index][3*i]));
743  // if(intersections_iterator->contains(p)){
744  // Ok, the target cell is in the intersection of the two regions
745  // Find source cell and get solution values
746  // 1. Extract the source cells close to the given point
747  std::vector<Mesh::IndexType> candidates;
748  GeoPrim::CBox little_box(source_boxes[source_region_index][1].around(p));
749  // weed out the cells by checking distance from centroid
750  global.FunctionEntry("WeedOut");
751  std::vector<Mesh::IndexType>::iterator srccellit = source_cells.begin();
752  while(srccellit != source_cells.end()){
753  GeoPrim::CPoint p2(&(source_cell_centers[source_region_index][3*(*srccellit-1)]));
754  if((p2-p).norm() <= source_edgelen)
755  candidates.push_back(*srccellit);
756  srccellit++;
757  }
758  global.FunctionExit("WeedOut");
759  if(candidates.empty()){
760  global.FunctionEntry("CollideCellsBox");
761  Mesh::CollideCellsWithBox(source_coords,sourcemesh,little_box,source_cells,candidates);
762  global.FunctionExit("CollideCellsBox");
763  }
764  else{
765  std::vector<Mesh::IndexType> candcopy(candidates);
766  candidates.resize(0);
767  global.FunctionEntry("CollideCellsBoxMini");
768  Mesh::CollideCellsWithBox(source_coords,sourcemesh,little_box,candcopy,candidates);
769  global.FunctionExit("CollideCellsBoxMini");
770  }
771  GeoPrim::CVector natc;
772  unsigned int source_cell_id = 0;
773  if(candidates.empty()){
774  // std::cout << "did not find candidates for target mesh point: (" << target_region_index+1
775  // << "," << i+1 << "," << p << ") with box: " << little_box << " in source mesh: "
776  // << source_region_index+1 << " with bounding box " << source_boxes[source_region_index][0]
777  // << std::endl;
778  global.FunctionEntry("FindPointInAllCells");
779  source_cell_id = Mesh::FindPointInCells(p,source_coords,sourcemesh,source_cells,natc);
780  global.FunctionExit("FindPointInAllCells");
781  }
782  else{
783  global.FunctionEntry("FindPointInCells");
784  source_cell_id = Mesh::FindPointInCells(p,source_coords,sourcemesh,candidates,natc);
785  global.FunctionExit("FindPointInCells");
786  }
787  if(source_cell_id > 0){
788  (*tsmit)[i].push_back(std::make_pair(source_region_index+1,source_cell_id));
789  intersection_cells_iterator->push_back(i+1);
790  cell_found[i] = true;
791  }
792  else{ // resort to something else terrible
793  // Ouf << "Target Point: " << p << " with box [" << little_box << "] " << std::endl
794  // << "Source Region: " << source_region_index+1 << " with box ["
795  // << source_boxes[source_region_index][0] << "]" << std::endl
796  // << "Source Elements: " << std::endl;
797  GeoPrim::CVector d(1000,1000,1000);
798  unsigned int closest_cell = 0;
799  std::vector<Mesh::IndexType> *thecells;
800  if(candidates.empty())
801  thecells = &source_cells;
802  else
803  thecells = &candidates;
804  // std::vector<Mesh::IndexType>::iterator ci = source_cells.begin();
805  // while(ci != source_cells.end())
806  std::vector<Mesh::IndexType>::iterator ci = thecells->begin();
807  while(ci != thecells->end())
808  {
809  source_cell_id = *ci++;
810  GeoPrim::CVector dist = p - GeoPrim::CPoint(&source_cell_centers[source_region_index][3*(source_cell_id-1)]);
811  if(dist.mag() < d.mag()){
812  d = dist;
813  closest_cell = source_cell_id;
814  }
815  }
816  // Ouf << "Closest cell (" << source_region_index+1 << "," << closest_cell << ") : ("
817  // << d << ")" << std::endl;
818  (*tsmit)[i].push_back(std::make_pair(source_region_index+1,closest_cell));
819  intersection_cells_iterator->push_back(i+1);
820  // if(candidates.empty()){
821  // }
822  // else {
823  // std::vector<Mesh::IndexType>::iterator ci = candidates.begin();
824  // while(ci != candidates.end())
825  // {
826  // unsigned int source_cell_id = *ci++;
827  // Ouf << "Element " << source_cell_id << ": " << std::endl;
828  // std::vector<Mesh::IndexType>::iterator ni = sourcemesh[source_cell_id-1].begin();
829  // while(ni != sourcemesh[source_cell_id-1].end())
830  // {
831  // unsigned int node_id = *ni++;
832  // Ouf << node_id << ": "
833  // << GeoPrim::C3Point(source_coords[node_id])
834  // << std::endl;
835  // }
836 
837  // }
838 
839  // }
840  }
841  }
842  }
843  intersection_cells_iterator++;
844  intersections_iterator++;
845  }
846 
847 
848  // std::vector<bool>::iterator cfi = cell_found.begin();
849  // unsigned int number_not_found = 0;
850  // while(cfi != cell_found.end())
851  // if(!*cfi++)
852  // number_not_found++;
853 
854  if(!((*target_boxes_iterator)[0] == intersection_box) && StdOut)
855  *StdOut << "Mesh Box: " << (*target_boxes_iterator)[0]
856  << " Intersecting region: " << intersection_box << std::endl;
857  tsmit++;
858  target_mesh_iterator++;
859  target_block_iterator++;
860  target_intersections_iterator++;
861  tint_cells_iterator++;
862  target_boxes_iterator++;
863  target_connectivity_iterator++;
864  }
865 
866  global.FunctionExit("FindIntersectionCells");
867  // }
868  comm.Barrier();
869  if(!rank && StdOut)
870  std::cout << std::endl << "Preparing to write solution." << std::endl;
871  global.FunctionEntry("Reporting");
872  unsigned int number_not_found = 0;
873  bool reporting = false;
874  tsmit = target_source_matches.begin();
875  while(tsmit != target_source_matches.end())
876  {
877  unsigned int target_region_index = tsmit - target_source_matches.begin();
878  std::vector<std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> > >::iterator tci = tsmit->begin();
879  while(tci != tsmit->end()){
880  unsigned int target_cell_index = tci - tsmit->begin();
881  if(tci->empty()){
882  number_not_found++;
883  if(debug && StdOut)
884  *StdOut << "(" << target_region_index+1 << "," << target_cell_index+1 << ") not found." << std::endl;
885  }
886  else {
887  if(debug && StdOut){
888  *StdOut << "(" << target_region_index+1 << "," << target_cell_index+1 << ") has "
889  << tci->size() << " matches." << std::endl
890  << "target point: "
891  << GeoPrim::C3Point(&target_region_cell_centers[target_region_index][3*target_cell_index]) << std::endl
892  << "source(s): " << std::endl;
893  }
894  if(tci->size() == 1)
895  target_source_cells[target_region_index][target_cell_index] = (*tci)[0];
896  else {
897  GeoPrim::CPoint p(&target_region_cell_centers[target_region_index][3*target_cell_index]);
898  std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> >::iterator tcit = tci->begin();
899  GeoPrim::CVector dist(1000,1000,1000);
900  std::pair<Mesh::IndexType,Mesh::IndexType> match;
901  while(tcit != tci->end()){
902  Mesh::IndexType source_region_index = tcit->first - 1;
903  Mesh::IndexType source_cell_index = tcit->second - 1;
904  GeoPrim::CVector d = p - GeoPrim::CPoint(&source_cell_centers[source_region_index][3*source_cell_index]);
905  if(d.mag() < dist.mag()){
906  dist = d;
907  match = *tcit;
908  }
909  tcit++;
910  }
911  target_source_cells[target_region_index][target_cell_index] = match;
912  }
913  std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> >::iterator tcit = tci->begin();
914  while(tcit != tci->end())
915  {
916  if(debug && StdOut){
917  *StdOut << "(" << tcit->first << "," << tcit->second << ") : ("
918  << GeoPrim::C3Point(&source_cell_centers[tcit->first-1][3*(tcit->second-1)])
919  << ") : ("
920  << GeoPrim::C3Point(&source_cell_centers[tcit->first-1][3*(tcit->second-1)]) -
921  GeoPrim::C3Point(&target_region_cell_centers[target_region_index][3*target_cell_index])
922  << ")" << std::endl;
923  }
924  tcit++;
925  }
926 
927  }
928  tci++;
929  }
930  tsmit++;
931  }
932  if(number_not_found > 0 && StdOut)
933  *StdOut << number_not_found << "/" << total_target_cells << " cells not found or processed." << std::endl;
934  global.FunctionExit("Reporting");
935  }
936  // Searching complete... Now do the solution transfer part
937 
938  // get rid of the rudding meshes
939  // sourcegrid.DestroyGrids();
940  // targetgrid.DestroyGrids();
941 
942  global.FunctionEntry("SolutionRead");
943  // read the source solution
944  if(sourcegrid.ReadAllSolutions(timestamp)){
945  std::cout << "Could not read source solutions." << std::endl;
946  exit(1);
947  }
948  global.FunctionExit("SolutionRead");
949 
950  // allocate space for the target solution
951  unsigned int ng = sourcegrid.Blocks()[0].NGhostLayers();
952  double unknown_number = sourcegrid.UnknownNumber();
953  double time = sourcegrid.Time();
954  targetgrid.SetGhostLayers(ng);
955  targetgrid.CreateSolutions();
956  std::vector<FloGridBlock> &targetgrids = targetgrid.Blocks();
957 
958  // Populate an array of flat source cell ids for performance
959  source_cell_ids.resize(NSourceRegions);
960  std::vector<FloGridBlock> &sourcegrids = sourcegrid.Blocks();
961  std::vector<FloGridBlock>::iterator sbi = sourcegrids.begin();
962  while(sbi != sourcegrids.end()){
963  unsigned int source_grid_index = sbi - sourcegrids.begin();
964  Mesh::BSExtent<Mesh::IndexType> real_extent;
965  Mesh::BSExtent<Mesh::IndexType> ghost_extent;
966  std::vector<Mesh::IndexType> source_block_extent(6,0);
967  std::vector<Mesh::IndexType> source_ghost_extent(6,0);
968  unsigned int isize = sbi->isize();
969  unsigned int jsize = sbi->jsize();
970  unsigned int ksize = sbi->ksize();
971  source_block_extent[0] = source_block_extent[2] = source_block_extent[4] = 1 + ng;
972  source_ghost_extent[0] = source_ghost_extent[2] = source_ghost_extent[4] = 1;
973  source_block_extent[1] = isize+ng;
974  source_block_extent[3] = jsize+ng;
975  source_block_extent[5] = ksize+ng;
976  source_ghost_extent[1] = isize+(2*ng);
977  source_ghost_extent[3] = jsize+(2*ng);
978  source_ghost_extent[5] = ksize+(2*ng);
979  real_extent.Init(source_block_extent);
980  ghost_extent.Init(source_ghost_extent);
981  ghost_extent.GetFlatIndices(real_extent,source_cell_ids[source_grid_index]);
982  sbi++;
983  }
984 
985  // Transfer the data
986  global.FunctionEntry("DataTransfer");
987  target_cell_ids.resize(NTargetRegions);
988  std::vector<std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> > >::iterator tsci = target_source_cells.begin();
989  while(tsci != target_source_cells.end()){
990  unsigned int target_region_index = tsci - target_source_cells.begin();
991  // std::vector<Mesh::IndexType> &target_cell_ids = all_target_cell_ids[target_region_index];
992  std::vector<Mesh::IndexType> &target_cell_ids_local = target_cell_ids[target_region_index];
993  Mesh::BSExtent<Mesh::IndexType> real_extent;
994  Mesh::BSExtent<Mesh::IndexType> ghost_extent;
995  std::vector<Mesh::IndexType> target_block_extent(6,0);
996  std::vector<Mesh::IndexType> target_ghost_extent(6,0);
997  unsigned int isize = targetgrids[target_region_index].isize();
998  unsigned int jsize = targetgrids[target_region_index].jsize();
999  unsigned int ksize = targetgrids[target_region_index].ksize();
1000  target_block_extent[0] = target_block_extent[2] = target_block_extent[4] = 1 + ng;
1001  target_ghost_extent[0] = target_ghost_extent[2] = target_ghost_extent[4] = 1;
1002  target_block_extent[1] = isize+ng;
1003  target_block_extent[3] = jsize+ng;
1004  target_block_extent[5] = ksize+ng;
1005  target_ghost_extent[1] = isize+(2*ng);
1006  target_ghost_extent[3] = jsize+(2*ng);
1007  target_ghost_extent[5] = ksize+(2*ng);
1008  real_extent.Init(target_block_extent);
1009  ghost_extent.Init(target_ghost_extent);
1010  // ghost_extent.GetFlatIndices(target_block_extent,target_cell_ids_local);
1011  ghost_extent.GetFlatIndices(real_extent,target_cell_ids_local);
1012  unsigned int target_block = (isize+(2*ng))*(jsize+(2*ng))*(ksize+(2*ng));
1013  // ghost_extent.GetFlatIndices(real_extent,all_target_cell_ids[target_region_index]);
1014  std::vector<std::pair<Mesh::IndexType,Mesh::IndexType> >::iterator tci = tsci->begin();
1015  while(tci != tsci->end()){
1016  unsigned int target_cell_index = target_cell_ids_local[tci - tsci->begin()] - 1;
1017  unsigned int source_region_index = tci->first - 1;
1018  unsigned int source_cell_index = source_cell_ids[source_region_index][tci->second -1] - 1;
1019  unsigned int source_block = (sourcegrids[source_region_index].isize()+(2*ng))*
1020  (sourcegrids[source_region_index].jsize()+(2*ng))*(sourcegrids[source_region_index].ksize()+(2*ng));
1021 
1022  targetgrids[target_region_index].Solution()[target_cell_index] =
1023  sourcegrids[source_region_index].Solution()[source_cell_index];
1024 
1025  targetgrids[target_region_index].Solution()[target_block+target_cell_index] =
1026  sourcegrids[source_region_index].Solution()[source_block+source_cell_index];
1027 
1028  targetgrids[target_region_index].Solution()[(2*target_block) + target_cell_index] =
1029  sourcegrids[source_region_index].Solution()[(2*source_block) + source_cell_index];
1030 
1031  targetgrids[target_region_index].Solution()[(3*target_block) + target_cell_index] =
1032  sourcegrids[source_region_index].Solution()[(3*source_block) + source_cell_index];
1033 
1034  targetgrids[target_region_index].Solution()[(4*target_block) + target_cell_index] =
1035  sourcegrids[source_region_index].Solution()[(4*source_block) + source_cell_index];
1036 
1037  tci++;
1038  }
1039  tsci++;
1040  }
1041  global.FunctionExit("DataTransfer");
1042  global.FunctionEntry("GhostZoneFlood");
1043  std::vector<FloGridBlock>::iterator tgi = targetgrids.begin();
1044  while(tgi != targetgrids.end()){
1045  unsigned int target_grid_index = tgi - targetgrids.begin();
1046  Mesh::BSExtent<Mesh::IndexType> full_extent;
1047  std::vector<Mesh::IndexType> flat_full_extent(6,0);
1048  unsigned int isize = tgi->isize();
1049  unsigned int jsize = tgi->jsize();
1050  unsigned int ksize = tgi->ksize();
1051  flat_full_extent[0] = flat_full_extent[2] = flat_full_extent[4] = 1;
1052  flat_full_extent[1] = isize+(2*ng);
1053  flat_full_extent[3] = jsize+(2*ng);
1054  flat_full_extent[5] = ksize+(2*ng);
1055  full_extent.Init(flat_full_extent);
1056  unsigned int target_block = (isize+(2*ng))*(jsize+(2*ng))*(ksize+(2*ng));
1057  bool do_iplane = true;
1058  if(do_iplane){
1059  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1060  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1061  flat_plane_extent[0] = flat_plane_extent[1] = ng+1;
1062  flat_plane_extent[2] = flat_plane_extent[4] = ng+1;
1063  flat_plane_extent[3] = jsize + ng;
1064  flat_plane_extent[5] = ksize + ng;
1065  plane_extent.Init(flat_plane_extent);
1066  std::vector<Mesh::IndexType> real_cell_ids;
1067  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1068  for(int ii = 1;ii <= ng;ii++){
1069  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1070  flat_ghost_layer[0] = flat_ghost_layer[1] = ii;
1071  flat_ghost_layer[2] = flat_plane_extent[2];
1072  flat_ghost_layer[3] = flat_plane_extent[3];
1073  flat_ghost_layer[4] = flat_plane_extent[4];
1074  flat_ghost_layer[5] = flat_plane_extent[5];
1075  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1076  std::vector<Mesh::IndexType> ghost_cell_ids;
1077  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1078  assert(ghost_cell_ids.size() == real_cell_ids.size());
1079  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1080  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1081  while(rci != real_cell_ids.end()){
1082  unsigned int real_cell_index = *rci++ - 1;
1083  unsigned int ghost_cell_index = *gci++ - 1;
1084  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1085  targetgrids[target_grid_index].Solution()[real_cell_index];
1086  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1087  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1088  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1089  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1090  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1091  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1092  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1093  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1094  }
1095  }
1096  }
1097  if(do_iplane){
1098  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1099  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1100  flat_plane_extent[0] = flat_plane_extent[1] = isize+ng;
1101  flat_plane_extent[2] = flat_plane_extent[4] = ng+1;
1102  flat_plane_extent[3] = jsize + ng;
1103  flat_plane_extent[5] = ksize + ng;
1104  plane_extent.Init(flat_plane_extent);
1105  std::vector<Mesh::IndexType> real_cell_ids;
1106  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1107  for(int ii = 1;ii <= ng;ii++){
1108  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1109  flat_ghost_layer[0] = flat_ghost_layer[1] = isize+ng+ii;
1110  flat_ghost_layer[2] = flat_plane_extent[2];
1111  flat_ghost_layer[3] = flat_plane_extent[3];
1112  flat_ghost_layer[4] = flat_plane_extent[4];
1113  flat_ghost_layer[5] = flat_plane_extent[5];
1114  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1115  std::vector<Mesh::IndexType> ghost_cell_ids;
1116  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1117  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1118  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1119  assert(ghost_cell_ids.size() == real_cell_ids.size());
1120  while(rci != real_cell_ids.end()){
1121  unsigned int real_cell_index = *rci++ - 1;
1122  unsigned int ghost_cell_index = *gci++ - 1;
1123  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1124  targetgrids[target_grid_index].Solution()[real_cell_index];
1125  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1126  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1127  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1128  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1129  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1130  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1131  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1132  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1133  }
1134  // while(rci != real_cell_ids.end())
1135  // *gci++ = *rci++;
1136  }
1137  }
1138  bool do_jplane = true;
1139  if(do_jplane){
1140  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1141  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1142  flat_plane_extent[0] = 1;
1143  flat_plane_extent[1] = isize+(2*ng);
1144  flat_plane_extent[2] = flat_plane_extent[3] = ng+1;
1145  flat_plane_extent[4] = ng+1;
1146  flat_plane_extent[5] = ksize + ng;
1147  plane_extent.Init(flat_plane_extent);
1148  std::vector<Mesh::IndexType> real_cell_ids;
1149  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1150  for(int ii = 1;ii <= ng;ii++){
1151  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1152  flat_ghost_layer[0] = flat_plane_extent[0];
1153  flat_ghost_layer[1] = flat_plane_extent[1];
1154  flat_ghost_layer[2] = flat_ghost_layer[3] = ii;
1155  flat_ghost_layer[4] = flat_plane_extent[4];
1156  flat_ghost_layer[5] = flat_plane_extent[5];
1157  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1158  std::vector<Mesh::IndexType> ghost_cell_ids;
1159  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1160  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1161  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1162  assert(ghost_cell_ids.size() == real_cell_ids.size());
1163  while(rci != real_cell_ids.end()){
1164  unsigned int real_cell_index = *rci++ - 1;
1165  unsigned int ghost_cell_index = *gci++ - 1;
1166  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1167  targetgrids[target_grid_index].Solution()[real_cell_index];
1168  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1169  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1170  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1171  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1172  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1173  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1174  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1175  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1176  }
1177  // while(rci != real_cell_ids.end())
1178  // *gci++ = *rci++;
1179  }
1180  }
1181  if(do_jplane){
1182  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1183  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1184  flat_plane_extent[0] = 1;
1185  flat_plane_extent[1] = isize+(2*ng);
1186  flat_plane_extent[2] = flat_plane_extent[3] = jsize+ng;
1187  flat_plane_extent[4] = ng+1;
1188  flat_plane_extent[5] = ksize + ng;
1189  plane_extent.Init(flat_plane_extent);
1190  std::vector<Mesh::IndexType> real_cell_ids;
1191  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1192  for(int ii = 1;ii <= ng;ii++){
1193  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1194  flat_ghost_layer[0] = flat_plane_extent[0];
1195  flat_ghost_layer[1] = flat_plane_extent[1];
1196  flat_ghost_layer[2] = flat_ghost_layer[3] = ii+jsize+ng;
1197  flat_ghost_layer[4] = flat_plane_extent[4];
1198  flat_ghost_layer[5] = flat_plane_extent[5];
1199  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1200  std::vector<Mesh::IndexType> ghost_cell_ids;
1201  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1202  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1203  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1204  assert(ghost_cell_ids.size() == real_cell_ids.size());
1205  while(rci != real_cell_ids.end()){
1206  unsigned int real_cell_index = *rci++ - 1;
1207  unsigned int ghost_cell_index = *gci++ - 1;
1208  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1209  targetgrids[target_grid_index].Solution()[real_cell_index];
1210  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1211  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1212  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1213  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1214  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1215  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1216  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1217  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1218  }
1219  // while(rci != real_cell_ids.end())
1220  // *gci++ = *rci++;
1221  }
1222  }
1223  bool do_kplane = true;
1224  if(do_kplane){
1225  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1226  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1227  flat_plane_extent[0] = 1;
1228  flat_plane_extent[1] = isize+(2*ng);
1229  flat_plane_extent[2] = 1;
1230  flat_plane_extent[3] = jsize+(2*ng);
1231  flat_plane_extent[4] = flat_plane_extent[5] = ng+1;
1232  plane_extent.Init(flat_plane_extent);
1233  std::vector<Mesh::IndexType> real_cell_ids;
1234  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1235  for(int ii = 1;ii <= ng;ii++){
1236  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1237  flat_ghost_layer[0] = flat_plane_extent[0];
1238  flat_ghost_layer[1] = flat_plane_extent[1];
1239  flat_ghost_layer[2] = flat_plane_extent[2];
1240  flat_ghost_layer[3] = flat_plane_extent[3];
1241  flat_ghost_layer[4] = flat_ghost_layer[5] = ii;
1242  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1243  std::vector<Mesh::IndexType> ghost_cell_ids;
1244  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1245  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1246  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1247  assert(ghost_cell_ids.size() == real_cell_ids.size());
1248  while(rci != real_cell_ids.end()){
1249  unsigned int real_cell_index = *rci++ - 1;
1250  unsigned int ghost_cell_index = *gci++ - 1;
1251  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1252  targetgrids[target_grid_index].Solution()[real_cell_index];
1253  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1254  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1255  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1256  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1257  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1258  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1259  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1260  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1261  }
1262  // while(rci != real_cell_ids.end())
1263  // *gci++ = *rci++;
1264  }
1265  }
1266  if(do_kplane){
1267  Mesh::BSExtent<Mesh::IndexType> plane_extent;
1268  std::vector<Mesh::IndexType> flat_plane_extent(6,0);
1269  flat_plane_extent[0] = 1;
1270  flat_plane_extent[1] = isize+(2*ng);
1271  flat_plane_extent[2] = 1;
1272  flat_plane_extent[3] = jsize+(2*ng);
1273  flat_plane_extent[4] = flat_plane_extent[5] = ng+ksize;
1274  plane_extent.Init(flat_plane_extent);
1275  std::vector<Mesh::IndexType> real_cell_ids;
1276  full_extent.GetFlatIndices(plane_extent,real_cell_ids);
1277  for(int ii = 1;ii <= ng;ii++){
1278  std::vector<Mesh::IndexType> flat_ghost_layer(6,0);
1279  flat_ghost_layer[0] = flat_plane_extent[0];
1280  flat_ghost_layer[1] = flat_plane_extent[1];
1281  flat_ghost_layer[2] = flat_plane_extent[2];
1282  flat_ghost_layer[3] = flat_plane_extent[3];
1283  flat_ghost_layer[4] = flat_ghost_layer[5] = ii+ng+ksize;
1284  Mesh::BSExtent<Mesh::IndexType> ghost_plane_extent(flat_ghost_layer);
1285  std::vector<Mesh::IndexType> ghost_cell_ids;
1286  full_extent.GetFlatIndices(ghost_plane_extent,ghost_cell_ids);
1287  std::vector<Mesh::IndexType>::iterator rci = real_cell_ids.begin();
1288  std::vector<Mesh::IndexType>::iterator gci = ghost_cell_ids.begin();
1289  assert(ghost_cell_ids.size() == real_cell_ids.size());
1290  while(rci != real_cell_ids.end()){
1291  unsigned int real_cell_index = *rci++ - 1;
1292  unsigned int ghost_cell_index = *gci++ - 1;
1293  targetgrids[target_grid_index].Solution()[ghost_cell_index] =
1294  targetgrids[target_grid_index].Solution()[real_cell_index];
1295  targetgrids[target_grid_index].Solution()[target_block+ghost_cell_index] =
1296  targetgrids[target_grid_index].Solution()[target_block+real_cell_index];
1297  targetgrids[target_grid_index].Solution()[(2*target_block)+ghost_cell_index] =
1298  targetgrids[target_grid_index].Solution()[(2*target_block)+real_cell_index];
1299  targetgrids[target_grid_index].Solution()[(3*target_block)+ghost_cell_index] =
1300  targetgrids[target_grid_index].Solution()[(3*target_block)+real_cell_index];
1301  targetgrids[target_grid_index].Solution()[(4*target_block)+ghost_cell_index] =
1302  targetgrids[target_grid_index].Solution()[(4*target_block)+real_cell_index];
1303  }
1304  // while(rci != real_cell_ids.end())
1305  // *gci++ = *rci++;
1306  }
1307  }
1308  tgi++;
1309  }
1310 
1311  global.FunctionExit("GhostZoneFlood");
1312 
1313  unsigned int ncells = 0;
1314  if(rank == 0) {
1315  std::vector<double> sourcemin(5,std::numeric_limits<double>::max());
1316  std::vector<double> sourcemax(5,-1.0*std::numeric_limits<double>::max());
1317  std::vector<double> sourcemean(5,0);
1318  std::vector<std::vector<Mesh::IndexType> >::iterator sri = source_cell_ids.begin();
1319  while(sri != source_cell_ids.end()){
1320  unsigned int source_region_index = sri - source_cell_ids.begin();
1321  std::vector<Mesh::IndexType>::iterator sci = sri->begin();
1322  unsigned int blocksize = (sourcegrids[source_region_index].isize()+(2*ng)) *
1323  (sourcegrids[source_region_index].jsize()+(2*ng)) * (sourcegrids[source_region_index].ksize()+(2*ng));
1324  while(sci != sri->end()){
1325  ncells++;
1326  unsigned int source_cell_index = *sci++-1;
1327 
1328  double val = sourcegrids[source_region_index].Solution()[source_cell_index];
1329  sourcemean[0] += val;
1330  if(val > sourcemax[0])
1331  sourcemax[0] = val;
1332  if(val < sourcemin[0])
1333  sourcemin[0] = val;
1334 
1335  val = sourcegrids[source_region_index].Solution()[blocksize+source_cell_index];
1336  sourcemean[1] += val;
1337  if(val > sourcemax[1])
1338  sourcemax[1] = val;
1339  if(val < sourcemin[1])
1340  sourcemin[1] = val;
1341 
1342  val = sourcegrids[source_region_index].Solution()[(2*blocksize)+source_cell_index];
1343  sourcemean[2] += val;
1344  if(val > sourcemax[2])
1345  sourcemax[2] = val;
1346  if(val < sourcemin[2])
1347  sourcemin[2] = val;
1348 
1349  val = sourcegrids[source_region_index].Solution()[(3*blocksize)+source_cell_index];
1350  sourcemean[3] += val;
1351  if(val > sourcemax[3])
1352  sourcemax[3] = val;
1353  if(val < sourcemin[3])
1354  sourcemin[3] = val;
1355 
1356  val = sourcegrids[source_region_index].Solution()[(4*blocksize)+source_cell_index];
1357  sourcemean[4] += val;
1358  if(val > sourcemax[4])
1359  sourcemax[4] = val;
1360  if(val < sourcemin[4])
1361  sourcemin[4] = val;
1362 
1363  }
1364  sri++;
1365  }
1366  if(ncells != total_source_cells)
1367  std::cout << "Warning: Source soln extents indicate only " << ncells
1368  << "/" << total_source_cells << " real cells." << std::endl;
1369  double scale = 1.0/static_cast<double>(ncells);
1370  sourcemean[0] *= scale;
1371  sourcemean[1] *= scale;
1372  sourcemean[2] *= scale;
1373  sourcemean[3] *= scale;
1374  sourcemean[4] *= scale;
1375  std::cout << "Source Data (min,max,mean):" << std::endl
1376  << "rho : (" << sourcemin[0] << "," << sourcemax[0]
1377  << "," << sourcemean[0] << ")" << std::endl
1378  << "rho-u: (" << sourcemin[1] << "," << sourcemax[1]
1379  << "," << sourcemean[1] << ")" << std::endl
1380  << "rho-v: (" << sourcemin[2] << "," << sourcemax[2]
1381  << "," << sourcemean[2] << ")" << std::endl
1382  << "rho-w: (" << sourcemin[3] << "," << sourcemax[3]
1383  << "," << sourcemean[3] << ")" << std::endl
1384  << "rho-E: (" << sourcemin[4] << "," << sourcemax[4]
1385  << "," << sourcemean[4] << ")" << std::endl;
1386  }
1387  std::vector<double> targetmin(5,std::numeric_limits<double>::max());
1388  std::vector<double> ttargetmin(5,std::numeric_limits<double>::max());
1389  std::vector<double> targetmax(5,-1*std::numeric_limits<double>::max());
1390  std::vector<double> ttargetmax(5,-1*std::numeric_limits<double>::max());
1391  std::vector<double> targetmean(5,0);
1392  std::vector<double> ttargetmean(5,0);
1393  ncells = 0;
1394  std::vector<std::vector<Mesh::IndexType> >::iterator sri = target_cell_ids.begin();
1395  while(sri != target_cell_ids.end()){
1396  unsigned int target_region_index = sri - target_cell_ids.begin();
1397  std::vector<Mesh::IndexType>::iterator sci = sri->begin();
1398  unsigned int blocksize = (targetgrids[target_region_index].isize()+(2*ng)) *
1399  (targetgrids[target_region_index].jsize()+(2*ng))*(targetgrids[target_region_index].ksize()+(2*ng));
1400  while(sci != sri->end()){
1401  ncells++;
1402  unsigned int target_cell_index = *sci++-1;
1403  double val = targetgrids[target_region_index].Solution()[target_cell_index];
1404  targetmean[0] += val;
1405  if(val > targetmax[0])
1406  targetmax[0] = val;
1407  if(val < targetmin[0])
1408  targetmin[0] = val;
1409  val = targetgrids[target_region_index].Solution()[blocksize+ target_cell_index];
1410  targetmean[1] += val;
1411  if(val > targetmax[1])
1412  targetmax[1] = val;
1413  if(val < targetmin[1])
1414  targetmin[1] = val;
1415  val = targetgrids[target_region_index].Solution()[2*blocksize+target_cell_index];
1416  targetmean[2] += val;
1417  if(val > targetmax[2])
1418  targetmax[2] = val;
1419  if(val < targetmin[2])
1420  targetmin[2] = val;
1421  val = targetgrids[target_region_index].Solution()[3*blocksize+target_cell_index];
1422  targetmean[3] += val;
1423  if(val > targetmax[3])
1424  targetmax[3] = val;
1425  if(val < targetmin[3])
1426  targetmin[3] = val;
1427  val = targetgrids[target_region_index].Solution()[4*blocksize +target_cell_index];
1428  targetmean[4] += val;
1429  if(val > targetmax[4])
1430  targetmax[4] = val;
1431  if(val < targetmin[4])
1432  targetmin[4] = val;
1433 
1434  }
1435  sri++;
1436  }
1437  if(ncells != total_target_cells)
1438  std::cout << "Warning: Target soln extents indicate only " << ncells
1439  << "/" << total_target_cells << " real cells." << std::endl;
1440 
1441  comm.Reduce(targetmean,ttargetmean,IRAD::Comm::DTDOUBLE,IRAD::Comm::SUMOP,0);
1442  comm.Reduce(targetmin,ttargetmin,IRAD::Comm::DTDOUBLE,IRAD::Comm::MINOP,0);
1443  comm.Reduce(targetmax,ttargetmax,IRAD::Comm::DTDOUBLE,IRAD::Comm::MAXOP,0);
1444  double scale = 1.0/static_cast<double>(ncells);
1445  targetmean[0] *= scale;
1446  targetmean[1] *= scale;
1447  targetmean[2] *= scale;
1448  targetmean[3] *= scale;
1449  targetmean[4] *= scale;
1450  scale = 1.0/static_cast<double>(totalcells);
1451  ttargetmean[0] *= scale;
1452  ttargetmean[1] *= scale;
1453  ttargetmean[2] *= scale;
1454  ttargetmean[3] *= scale;
1455  ttargetmean[4] *= scale;
1456  if(verblevel > 0 && StdOut){
1457  if( (nproc==1) || ((rank > 0) && debug)){
1458  *StdOut << "Target Data (min,max,mean):" << std::endl
1459  << "rho : (" << targetmin[0] << "," << targetmax[0]
1460  << "," << targetmean[0] << ")" << std::endl
1461  << "rho-u: (" << targetmin[1] << "," << targetmax[1]
1462  << "," << targetmean[1] << ")" << std::endl
1463  << "rho-v: (" << targetmin[2] << "," << targetmax[2]
1464  << "," << targetmean[2] << ")" << std::endl
1465  << "rho-w: (" << targetmin[3] << "," << targetmax[3]
1466  << "," << targetmean[3] << ")" << std::endl
1467  << "rho-E: (" << targetmin[4] << "," << targetmax[4]
1468  << "," << targetmean[4] << ")" << std::endl;
1469  }
1470  else {
1471  if(debug){
1472  *StdOut << "Target Data (min,max,mean):" << std::endl
1473  << "rho : (" << targetmin[0] << "," << targetmax[0]
1474  << "," << targetmean[0] << ")" << std::endl
1475  << "rho-u: (" << targetmin[1] << "," << targetmax[1]
1476  << "," << targetmean[1] << ")" << std::endl
1477  << "rho-v: (" << targetmin[2] << "," << targetmax[2]
1478  << "," << targetmean[2] << ")" << std::endl
1479  << "rho-w: (" << targetmin[3] << "," << targetmax[3]
1480  << "," << targetmean[3] << ")" << std::endl
1481  << "rho-E: (" << targetmin[4] << "," << targetmax[4]
1482  << "," << targetmean[4] << ")" << std::endl;
1483  }
1484 
1485  *StdOut << "Target Data (min,max,mean):" << std::endl
1486  << "rho : (" << ttargetmin[0] << "," << ttargetmax[0]
1487  << "," << ttargetmean[0] << ")" << std::endl
1488  << "rho-u: (" << ttargetmin[1] << "," << ttargetmax[1]
1489  << "," << ttargetmean[1] << ")" << std::endl
1490  << "rho-v: (" << ttargetmin[2] << "," << ttargetmax[2]
1491  << "," << ttargetmean[2] << ")" << std::endl
1492  << "rho-w: (" << ttargetmin[3] << "," << ttargetmax[3]
1493  << "," << ttargetmean[3] << ")" << std::endl
1494  << "rho-E: (" << ttargetmin[4] << "," << ttargetmax[4]
1495  << "," << ttargetmean[4] << ")" << std::endl;
1496  }
1497 
1498  }
1499 
1500  global.FunctionEntry("SolutionWrite");
1501  if(rank == 0)
1502  targetgrid.OpenSolutionFile(timestamp,time,unknown_number);
1503  std::ostringstream Ostr;
1504  Ostr << std::scientific << std::setprecision(16);
1505  int datasize = targetgrid.WriteBlocks(Ostr);
1506  std::vector<int> datasizes(nproc,0);
1507  comm.Gather(datasize,datasizes);
1508  if(rank > 0)
1509  comm._Send(const_cast<char *>(Ostr.str().c_str()),datasize,0,0);
1510  if(!rank){
1511  std::ofstream &Ouf = targetgrid.SolnFile();
1512  targetgrid.WriteBlocks(Ouf);
1513  for(unsigned int npi = 1;npi < nproc;npi++)
1514  {
1515  std::vector<char> recvbuf(datasizes[npi]);
1516  comm._Recv(&recvbuf[0],datasizes[npi],npi,0);
1517  std::vector<char>::iterator rbi = recvbuf.begin();
1518  while(rbi != recvbuf.end())
1519  Ouf << *rbi++;
1520  }
1521  targetgrid.CloseSolutionFile();
1522  }
1523  comm.Barrier();
1524  // targetgrid.WriteSolutionFile(timestamp,time,unknown_number);
1525  global.FunctionExit("SolutionWrite");
1526 
1527 // tgi = targetgrids.begin();
1528 // while(tgi != targetgrids.end()){
1529 // Mesh::NodalCoordinates nc;
1530 // unsigned int nnodes = tgi->NNodes();
1531 // nc.init(nnodes,&((tgi->Coords())[0]));
1532 // unsigned int block_index = tgi - targetgrids.begin();
1533 // unsigned int blocksize = (targetgrids[block_index].isize()+(2*ng)) *
1534 // (targetgrids[block_index].jsize()+(2*ng))*(targetgrids[block_index].ksize()+(2*ng));
1535 // std::ostringstream Ostr;
1536 // Ostr << "target_" << block_index+1 << ".vtk";
1537 // std::vector<double> emptyv;
1538 // emptyv.resize(0);
1539 // writeVtkData(nc,target_mesh[block_index],Ostr.str(),tgi->Solution(),target_cell_ids[block_index],blocksize);
1540 // tgi++;
1541 // }
1542 // tgi = sourcegrids.begin();
1543 // while(tgi != sourcegrids.end()){
1544 // Mesh::NodalCoordinates nc;
1545 // unsigned int nnodes = tgi->NNodes();
1546 // nc.init(nnodes,&((tgi->Coords())[0]));
1547 // unsigned int block_index = tgi - sourcegrids.begin();
1548 // unsigned int blocksize = (sourcegrids[block_index].isize()+(2*ng)) *
1549 // (sourcegrids[block_index].jsize()+(2*ng))*(sourcegrids[block_index].ksize()+(2*ng));
1550 // std::ostringstream Ostr;
1551 // Ostr << "source_" << block_index+1 << ".vtk";
1552 // std::vector<double> emptyv;
1553 // emptyv.resize(0);
1554 // writeVtkData(nc,source_mesh[block_index],Ostr.str(),tgi->Solution(),source_cell_ids[block_index],blocksize);
1555 // tgi++;
1556 // }
1557 
1558  // test();
1559  // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1560  // This stuff goes at the very end of the program, don't modify
1561  // below here.
1562  // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1563 
1564  global.FunctionExit("main");
1565  comm.Barrier();
1566  if(LogOut && debug)
1567  *LogOut << "All processors made it to the end." << std::endl;
1568  global.Finalize();
1569  // global.Profiler.Finalize();
1570  if(StdOut)
1571  // global.Report(Profiler.SummarizeSerialExecution(*StdOut);
1572  global.Report(*StdOut);
1573 
1574  return(0);
1575 }
Mesh::GetMeshCentroids
int GetMeshCentroids(Mesh::NodalCoordinates &nc, Mesh::Connectivity &conn, std::vector< double > &centroids)
Mesh::Connectivity
General connectivity object.
Definition: Mesh.H:334
d
const NT & d
Definition: rational_rotation.h:73
GeoPrim::CBox::P2
const CPoint & P2() const
Definition: GeoPrimitives.H:1184
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
intersection
CGAL_END_NAMESPACE CGAL_BEGIN_NAMESPACE Object intersection(const Line_2< R > &line1, const Line_2< R > &line2)
Definition: Line_2_Line_2_intersection.h:140
nvc::norm
T norm(const NVec< DIM, T > &v)
Definition: Rocon/include/NVec.h:195
Mesh::GetMeshBoxes
void GetMeshBoxes(const NodalCoordinates &nc, const Connectivity &ec, GeoPrim::CBox &mesh_box, GeoPrim::CBox &small_box, GeoPrim::CBox &large_box)
Bounding boxes for a mesh.
GeoPrim::CVector::mag
double mag() const
Definition: GeoPrimitives.H:289
Mesh::BSExtent::GetFlatIndices
void GetFlatIndices(const BSExtent< T > &extent, std::vector< T > &indices)
Definition: BSMesh.H:139
i
blockLoc i
Definition: read.cpp:79
Mesh::CollideCellsWithBox
int CollideCellsWithBox(Mesh::NodalCoordinates &nc, Mesh::Connectivity &conn, GeoPrim::CBox &box, std::vector< Mesh::IndexType > &candidates, std::vector< Mesh::IndexType > &cells)
Mesh::FindPointInCells
Mesh::IndexType FindPointInCells(const GeoPrim::CPoint &p, const NodalCoordinates &nc, const Connectivity &ec, const std::vector< Mesh::IndexType > &elements, GeoPrim::CVector &natc)
Locate element containing given physical point.
cimg_library::cimg::time
unsigned long time()
Get the value of a system timer with a millisecond precision.
Definition: CImg.h:4605
Mesh::BSExtent::Init
void Init(const std::vector< T > &inflatextent)
Definition: BSMesh.H:75
GeoPrim::CBox::merge
void merge(const CBox &b)
Definition: GeoPrimitives.H:1291
Mesh::CollideMeshWithBox
int CollideMeshWithBox(Mesh::NodalCoordinates &nc, Mesh::Connectivity &conn, GeoPrim::CBox &box, std::vector< Mesh::IndexType > &cells)
rank
static int rank
Definition: advectest.C:66
dist
long double dist(long double *coord1, long double *coord2, int size)
Definition: utilities/rc_cmp/util.cpp:305
GeoPrim::C3Point
C3Vector C3Point
Definition: GeoPrimitives.H:763
Mesh::IndexType
IRAD::Primitive::IndexType IndexType
Definition: Mesh.H:57
Mesh::BSExtent
Simple Block Structured Mesh object.
Definition: BSMesh.H:16
GeoPrim::CBox::empty
bool empty() const
Definition: GeoPrimitives.H:1185
GeoPrim::CBox::P1
const CPoint & P1() const
Definition: GeoPrimitives.H:1183

Here is the call graph for this function:

void test ( void  )

Definition at line 99 of file flotsam.C.

References Mesh::FindPointInCells(), CPoint::init(), Mesh::NodalCoordinates::init_node(), and Mesh::Connectivity::Sync().

Referenced by node::intersect(), GridConversion::ParallelTest(), Rocstar::ParallelTest(), Rocstar::Test(), GridConversion::Test(), and TRAIL_RemeshWrite().

100 {
101  Mesh::Connectivity conn;
102  Mesh::NodalCoordinates nc(24);
103  GeoPrim::CPoint p;
104  nc.init_node(1,p.init(.58992,.107061,0));
105  nc.init_node(2,p.init(.5946,0.103897,0));
106  nc.init_node(3,p.init(.5946,0.107263,0));
107  nc.init_node(4,p.init(.58992,0.110554,0));
108  nc.init_node(5,p.init(.58992,0.107061,.0508));
109  nc.init_node(6,p.init(.5946,.103897,0.0508));
110  nc.init_node(7,p.init(.5946,.107263,0.0508));
111  nc.init_node(8,p.init(.58992,0.110554,0.0508));
112  nc.init_node(9,p.init(0.5946,0.100531,0));
113  nc.init_node(10,p.init(0.59928,0.0974934,0));
114  nc.init_node(11,p.init(0.59928,0.100733 ,0));
115  nc.init_node(12,p.init(0.5946,0.103897,0));
116  nc.init_node(13,p.init(0.5946,0.100531,0.0508));
117  nc.init_node(14,p.init(0.59928,0.0974934,0.0508));
118  nc.init_node(15,p.init(0.59928,0.100733,0.0508));
119  nc.init_node(16,p.init(0.5946,0.103897,0.0508));
120  nc.init_node(17,p.init(0.699524,-0.0237386,0));
121  nc.init_node(18,p.init(0.694296,-0.0245241,0));
122  nc.init_node(19,p.init(0.694257,-0.0253172,0));
123  nc.init_node(20,p.init(0.699457,-0.0244292,0));
124  nc.init_node(21,p.init(0.699524,-0.0237386,0.0508));
125  nc.init_node(22,p.init(0.694296,-0.0245241,0.0508));
126  nc.init_node(23,p.init(0.694257,-0.0253172,0.0508));
127  nc.init_node(24,p.init(0.699457,-0.0244292,0.0508));
128  conn.resize(3);
129  conn[0].push_back(1);
130  conn[0].push_back(2);
131  conn[0].push_back(3);
132  conn[0].push_back(4);
133  conn[0].push_back(5);
134  conn[0].push_back(6);
135  conn[0].push_back(7);
136  conn[0].push_back(8);
137  conn[1].push_back(9);
138  conn[1].push_back(10);
139  conn[1].push_back(11);
140  conn[1].push_back(12);
141  conn[1].push_back(13);
142  conn[1].push_back(14);
143  conn[1].push_back(15);
144  conn[1].push_back(16);
145  conn[2].push_back(17);
146  conn[2].push_back(18);
147  conn[2].push_back(19);
148  conn[2].push_back(20);
149  conn[2].push_back(21);
150  conn[2].push_back(22);
151  conn[2].push_back(23);
152  conn[2].push_back(24);
153  conn.Sync();
154 // Point: 0.595836 0.111593 0.0254
155 // Point: 0.696963 -0.0204658 0.0762
156  std::vector<Mesh::IndexType> els(3);
157  els[0] = 1;
158  els[1] = 2;
159  els[2] = 3;
160  GeoPrim::CVector natc;
161  unsigned int cell_id = 0;
162  std::vector<GeoPrim::CPoint> test_points(3);
163  test_points[0].init(.595836,.111593,0.0254);
164  test_points[1].init( 0.60162,0.103916,0.0254 );
165  test_points[2].init(0.696963 ,-0.0204658 ,0.0762);
166  std::vector<GeoPrim::CPoint>::iterator tpi = test_points.begin();
167  while(tpi != test_points.end())
168  cell_id = Mesh::FindPointInCells(*tpi++,nc,conn,els,natc);
169 
170 }
Mesh::Connectivity::Sync
void Sync()
Definition: Mesh.C:246
GeoPrim::CPoint::init
CPoint & init()
Definition: GeoPrimitives.H:52
Mesh::Connectivity
General connectivity object.
Definition: Mesh.H:334
Mesh::FindPointInCells
Mesh::IndexType FindPointInCells(const GeoPrim::CPoint &p, const NodalCoordinates &nc, const Connectivity &ec, const std::vector< Mesh::IndexType > &elements, GeoPrim::CVector &natc)
Locate element containing given physical point.

Here is the call graph for this function:

Here is the caller graph for this function:

void test2 ( void  )

Definition at line 78 of file flotsam.C.

References Mesh::BSExtent< T >::GetFlatIndices(), and Mesh::BSExtent< T >::Init().

79 {
80  Mesh::BSExtent<Mesh::IndexType> ghost_extent;
81  Mesh::BSExtent<Mesh::IndexType> real_extent;
82  std::vector<Mesh::IndexType> flat_ghost(6,0);
83  std::vector<Mesh::IndexType> flat_real(6,0);
84  flat_ghost[0] = flat_ghost[2] = flat_ghost[4] = 1;
85  flat_ghost[1] = flat_ghost[3] = flat_ghost[5] = 4;
86  flat_real[0] = flat_real[2] = flat_real[4] = 2;
87  flat_real[1] = flat_real[3] = flat_real[5] = 3;
88  ghost_extent.Init(flat_ghost);
89  real_extent.Init(flat_real);
90  std::vector<Mesh::IndexType> flat_real_indices;
91  ghost_extent.GetFlatIndices(real_extent,flat_real_indices);
92  std::vector<Mesh::IndexType>::iterator fri = flat_real_indices.begin();
93  while(fri != flat_real_indices.end()){
94  std::cout << *fri++ << " ";
95  }
96  std::cout << std::endl;
97 }
Mesh::BSExtent::GetFlatIndices
void GetFlatIndices(const BSExtent< T > &extent, std::vector< T > &indices)
Definition: BSMesh.H:139
Mesh::BSExtent::Init
void Init(const std::vector< T > &inflatextent)
Definition: BSMesh.H:75
Mesh::BSExtent
Simple Block Structured Mesh object.
Definition: BSMesh.H:16

Here is the call graph for this function:

void writeVtkData ( Mesh::NodalCoordinates &  nc,
Mesh::Connectivity &  con,
const std::string &  filename,
std::vector< double > &  soln,
std::vector< Mesh::IndexType >  indices,
unsigned int  stride 
)

Definition at line 172 of file flotsam.C.

References Mesh::Connectivity::Esize(), i, j, Mesh::NodalCoordinates::NNodes(), and Mesh::Connectivity::SyncSizes().

176 {
177 
178  // std::cout << "Writing VTK data...";
179  // std::cout.flush( );
180 
181  std::ofstream ofs;
182  ofs.open( filename.c_str());
183 
184  if( !ofs.is_open( ) )
185  {
186  std::cerr << "Cannot write VTK file: " << filename << std::endl;
187  std::cerr << "File: " << __FILE__ << std::endl;
188  std::cerr << "Line: " << __LINE__ << std::endl;
189  assert( false );
190  }
191 
192  unsigned int nnodes = nc.NNodes();
193  /* STEP 1: Write the header */
194  ofs << "# vtk DataFile Version 3.0\n";
195  ofs << "Global Mesh Output" << std::endl;
196  ofs << "ASCII\n";
197  ofs << "DATASET UNSTRUCTURED_GRID\n";
198  ofs << "POINTS " << nnodes << " double\n";
199 
200  /* STEP 2: Write the nodes */
201  for( int i=1; i <= nnodes; ++i )
202  {
203  ofs << GeoPrim::C3Point(nc[i]) << std::endl;
204  }
205 
206  /* STEP 3: Write the max element size */
207  int maxnodes = 0; /* the maximum number of nodes per element */
208  unsigned int nelem = con.size();
209  con.SyncSizes();
210  for( int i=0; i < nelem; ++i )
211  {
212  unsigned int esize = con[i].size();
213  if(maxnodes < esize)
214  maxnodes = esize;
215  }
216 
217  /* STEP 4: Write the element connectivity */
218  ofs << "CELLS " << nelem << " " << nelem*( maxnodes+1) << "\n";
219  for( int i=0; i < nelem; ++i )
220  {
221  ofs << con.Esize(i+1) << " ";
222  for( int j=0; j < con.Esize(i+1); ++j )
223  ofs << con[ i ][ j ]-1 << " ";
224  ofs << std::endl;
225  }
226 
227  /* STEP 5: Write the cell types */
228  ofs << "CELL_TYPES " << nelem << std::endl;
229  for( int i=0; i < nelem; ++i )
230  {
231 
232  if( con.Esize(i+1) == 8 )
233  ofs << "12\n";
234  else if( con.Esize(i+1) == 4 )
235  ofs << "10\n";
236  else {
237  std::cerr << "Undefined element type!\n";
238  std::cerr << "File: " << __FILE__ << std::endl;
239  std::cerr << "Line: " << __LINE__ << std::endl;
240  }
241 
242  }
243 
244  /* STEP 6: Write the solution data attached to the points */
245  bool wroteHeader = false;
246  std::vector<std::string> Name(5);
247  Name[0] = "rho";
248  Name[1] = "rho-u";
249  Name[2] = "rho-v";
250  Name[3] = "rho-w";
251  Name[4] = "rho-E";
252  if(!soln.empty()){
253  for( int i=0; i < 5; ++i )
254  {
255  if( !wroteHeader )
256  {
257  ofs << "CELL_DATA " << nelem << std::endl;
258  wroteHeader = true;
259  }
260  ofs << "SCALARS " << Name[i] << " double" << std::endl;
261  ofs << "LOOKUP_TABLE default\n";
262  for( int j=0; j < nelem; ++j )
263  {
264 
265  ofs << soln[ i*stride + indices[j] - 1] << " ";
266  }
267  ofs << std::endl;
268  }
269  }
270  ofs.close( );
271  // std::cout << "[DONE]\n";
272 }
Mesh::Connectivity::Esize
IndexType Esize(IndexType n) const
Definition: Mesh.H:365
i
blockLoc i
Definition: read.cpp:79
Mesh::NodalCoordinates::NNodes
int NNodes() const
Definition: Mesh.H:278
Mesh::Connectivity::SyncSizes
void SyncSizes()
Definition: Mesh.C:281
j
j indices j
Definition: Indexing.h:6
GeoPrim::C3Point
C3Vector C3Point
Definition: GeoPrimitives.H:763

Here is the call graph for this function: