src/expde/grid/pridxpa.cc

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
//
// print_pa.cc
//
// ------------------------------------------------------------



// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <time.h>
 #include <math.h>
 #define Problemzeile  Datei->setf(ios::scientific,ios::floatfield);       
#else
 #include <iostream>
 #include <fstream>
 #include <ctime>
 #include <cmath>
 #define Problemzeile  Datei->setf(std::ios::scientific,std::ios::floatfield);
#endif 
  
// Include level 0:
#include "../parser.h"

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"
 
// Include level 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"

// Include level 5:
#include "gpar.h"
#include "parallel.h"
#include "mgcoeff.h"
#include "sto_man.h"
#include "gridbase.h"
#include "grid.h"


// 1. coarese level
// 2. finest  level
// 3. region  level
// 4. surface level
// 5. cells




// 2. finest level





void Grid_base::Print_Variable_OpenDx_parallel(ofstream *Datei, 
                                               int number_var) {
  int number_cells, number, ebene;
  int number_points;
  int n_proc,i;

  int *number_cells_process;
  int *number_points_process;
  int total_number_cells;
  int total_number_points;

  MPI_Status status;

  D3vector loc;
  int num, start_num, start_poi;

  int max_number_cells;
  int max_number_points;

  int*    buffer_send_int;
  double* buffer_send_double;

  int*    buffer_rec_int;
  double* buffer_rec_double;

  Index3D ind;
  int number_of_vars;   // wird in der Kopfzeile eingesetzt

  number_of_vars=1;

  MPI_Barrier(comm);

  if(Is_storage_initialized()==false) { 
    cout << " \n Fehler in Print_Variable_OpenDx! Initialisierung fehlt! "
         << endl;
    return;
  }
  if(Give_max_num_var() <= number_var) {
    cout << " \n Fehler in Print_Variable_OPENDX! number_var zu gross! " 
         << endl;
  }
  else { 
    // Teil 0: Write information
    *Datei << "# File created by ExPDE " << endl;
    *Datei << "# dx file format for OpenDx " << endl;
    *Datei << "# scalar value " << endl;

    
    // Teil 1: Anzahl der Zellen zaehlen
    // Teil 1.1: Innerer Zellen
    ind = my_index;
    number_cells = Recursion_Count_Cells(ind);
    // Teil 1.2: Randzellen
    iterate_hash2 {
      if(bocell->Give_Index()<<my_index) 
        number_cells = opendx_bo_cell(bocell,Datei,false,number_cells);
    }
    // Teil 1.3: Anzahl der Punkte
    // a) innen und randnah
    number_points=0;
    iterate_hash1 {

      /*
if(point1->ind == Index3D(9,9,2)) {
  cout << " Hi: sort: " 
       << " \n  my rank: " << my_rank 
       << " \n  typ: " << point1->typ
       << " \n  responsible? " << (point1->Give_Index()<<my_index)
       << endl;
  point1->ind.Print();
  cout << endl;
  my_index.Print();
  cout << endl;
}
      */


      if((point1->typ >= interior || point1->typ==parallel_p) && 
         point1->Give_Index()<<my_index && 
         point1->finest_level == max_level) { 
        point1->nummer = number_points; 
        number_points++; 
      }
      else point1->nummer = -1;
    }
    // b) Randpunkte
    iterate_hash3 {
      bo2point->Set_number_avs(number_points);
      ++number_points;
    }
    // c) Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Give_Index()<<my_index) {
        if(bocell->Exists_bocellpoint()) {
          bocell->Set_number_avs(number_points);
          number_points++;
        }
      }
    }

    // Teil 1.5: gather data 
    n_proc = give_number_of_processes();  
    number_cells_process = new int[n_proc];  
    number_points_process = new int[n_proc];  
    MPI_Gather(&number_points,1,MPI_INT,number_points_process,1,MPI_INT,  
               0,comm);  
    MPI_Gather(&number_cells,1,MPI_INT,number_cells_process,1,MPI_INT,  
               0,comm);  
    total_number_points =0;  
    total_number_cells =0;  
    max_number_points =0;  
    max_number_cells =0;  
    for(i=0;i<n_proc;++i) {  
      total_number_points = total_number_points + number_points_process[i];  
      total_number_cells  = total_number_cells  + number_cells_process[i];  
    }  
    for(i=0;i<n_proc;++i) {  
      if(number_points_process[i] > max_number_points)  
        max_number_points = number_points_process[i];  
      if(number_cells_process[i] > max_number_cells)  
        max_number_cells = number_cells_process[i];  
    }  
    // Teil 1.7: make buffer
    if(my_rank==0) {  
      buffer_rec_int    = new int[max_number_cells  * 4];  
      buffer_rec_double = new double[max_number_points * 3];  
    }  
    buffer_send_int    = new int[number_cells  * 4];  
    buffer_send_double = new double[number_points * 3];  

    // Teil 1.8: make buffer 
    if(my_rank==0) {  
      *Datei << "object 1 class array type float rank 1 shape 3 items "
             << total_number_points << " data follows" << endl;
    }


    // Teil 2: Koordinaten der Punkte eingeben 
    // a) set values in buffer 
    if(I_am_active()) {  
      // i) alle inneren Punkte  
      iterate_hash1 {  
        if(point1->nummer != -1) { 
          loc = transform_coord(point1->ind.coordinate());  
          buffer_send_double[point1->nummer*3]   = loc.x;  
          buffer_send_double[point1->nummer*3+1] = loc.y;  
          buffer_send_double[point1->nummer*3+2] = loc.z;  
        }  
      }  
      // ii) alle Randpunkte 
      iterate_hash3 {  
        loc = bo2point->transform_coord(Give_A(),H_mesh());  
        buffer_send_double[bo2point->nummer*3]   = loc.x;  
        buffer_send_double[bo2point->nummer*3+1] = loc.y;  
        buffer_send_double[bo2point->nummer*3+2] = loc.z;  
      }  
      // iii) alle Randzellfreiheitsgrad  
      iterate_hash2 {  
        if(bocell->Give_Index()<<my_index) {  
          if(bocell->Exists_bocellpoint()) {  
            loc =  transform_coord(bocell->ind.coordinate()) +   
              bocell->Local_coord_bocellpoint();  
            buffer_send_double[bocell->Number_avs()*3]   = loc.x;  
            buffer_send_double[bocell->Number_avs()*3+1] = loc.y;  
            buffer_send_double[bocell->Number_avs()*3+2] = loc.z;  
          }  
        }  
      }  
    }  
    // b) send 
    if(my_rank!=0) {  
      MPI_Send(buffer_send_double,number_points * 3,  
               MPI_DOUBLE,0,50,comm);  
    }  
    else {  
      start_num=0;  
      for(i=0;i<n_proc;++i) {  
        // c) receive  
        if(i!=0) {  
          MPI_Recv(buffer_rec_double,number_points_process[i] * 3,  
                   MPI_DOUBLE,i,50,comm,&status);  
        }  
        else {  
          for(num=0;num<number_points_process[i];++num) {  
            buffer_rec_double[num*3]   = buffer_send_double[num*3];  
            buffer_rec_double[num*3+1] = buffer_send_double[num*3+1];  
            buffer_rec_double[num*3+2] = buffer_send_double[num*3+2];  
          }  
        }  
        // d) print in file  
        for(num=0;num<number_points_process[i];++num) {  
          loc.x =  buffer_rec_double[num*3];  
          loc.y =  buffer_rec_double[num*3+1];  
          loc.z =  buffer_rec_double[num*3+2];  
          *Datei << " " << loc.x    
                 << " " << loc.y    
                 << " " << loc.z  
                 << "\n";  
        }  
        start_num = start_num + number_points_process[i];  
      }  
    }  


    // Teil 3: Zellen ausgeben 
    // Teil 3.0: Header file for tetrahedra:
    if(my_rank==0) {  
      *Datei << "object 2 class array type int rank 1 shape 4 items "
             << total_number_cells << " data follows" << endl;
    }

    // a) set values in buffer 
    if(I_am_active()) { 
      // Teil 3.1: Innere Zellen ausgeben 
      ind = my_index;  
      number=Recursion_Cells_AVS_parallel(ind,0,buffer_send_int);  
      // Teil 3.2: Randzellen 
      iterate_hash2 {  
        if(bocell->Give_Index()<<my_index) {
          number = avs_bo_cell_parallel(bocell,number,buffer_send_int);  
        }
      }
    }
    // b) send 
    if(my_rank!=0) {   
      MPI_Send(buffer_send_int,number_cells * 4,  
               MPI_INT,0,51,comm);  
    }  
    else {  
      start_num=0;  
      start_poi=0;  
      for(i=0;i<n_proc;++i) {  
        // c) receive   
        if(i!=0) {  
          MPI_Recv(buffer_rec_int,number_cells_process[i] * 4,  
                   MPI_INT,i,51,comm,&status);  
        }  
        else {  
          for(num=0;num<number_cells_process[i];++num) {  
            buffer_rec_int[num*4]   = buffer_send_int[num*4];   
            buffer_rec_int[num*4+1] = buffer_send_int[num*4+1];   
            buffer_rec_int[num*4+2] = buffer_send_int[num*4+2];   
            buffer_rec_int[num*4+3] = buffer_send_int[num*4+3];   
          }  
        }  
        // d) print in file  
        if(i<n_proc) {  
          for(num=0;num<number_cells_process[i];++num) {  
            *Datei << start_poi + buffer_rec_int[num*4]   
                   << " " << start_poi + buffer_rec_int[num*4+1]   
                   << " " << start_poi + buffer_rec_int[num*4+2]   
                   << " " << start_poi + buffer_rec_int[num*4+3]   
                   << "\n";  
          }  
        }  
        start_num = start_num + number_cells_process[i];  
        start_poi = start_poi + number_points_process[i];  
      }  
    }  
    if(my_rank==0) {
      *Datei << "attribute " << '"' << "element type" 
             << '"' << " string " << '"' << "tetrahedra" << '"' << "" << endl;
      *Datei << "attribute " << '"' << "ref" << '"' 
             << " string " << '"' << "positions" << '"' << "" << endl;
      
      // Teil 4: Header file for one variable:
      *Datei << "object 3 class array type float rank 0 items "
             << total_number_points << " data follows" << endl;
    }

    if(I_am_active()) {
    // Teil 2: Werte der Punkte ausgeben
    // a) set values in buffer
    // i) alle inneren Punkte
    iterate_hash1 {
      if(point1->nummer != -1) {
        ebene = Give_finest_level(point1->ind);
        buffer_send_double[point1->nummer] = 
          Give_variable(point1->ind,ebene)[number_var];


        /*      
if(Give_variable(point1->ind,ebene)[number_var] < 0.5) {
  if(point1->ind == Index3D(2,4,2)) {
    cout << " Hi: dx: " 
         <<  Give_variable(point1->ind,ebene)[number_var] 
         << " my rank: " << my_rank 
         << endl;
    point1->ind.Print();
    cout << endl;

    buffer_send_double[point1->nummer] = 666;
  }
}
        */
        



      }
    }
    // ii) alle Randpunkte
    iterate_hash3 {
      buffer_send_double[bo2point->nummer] = 
        Give_variable(bo2point->ind,bo2point->direction)[number_var];
    }
    // iii alle Randzellfreiheitsgrad
    iterate_hash2 {
      if(bocell->Give_Index()<<my_index) {
        if(bocell->Exists_bocellpoint()) {
          buffer_send_double[bocell->Number_avs()] =
            bocell->Give_var()[number_var];
        }
      }
    }
    }

    // b) send
    if(my_rank!=0) {
      MPI_Send(buffer_send_double,number_points,
               MPI_DOUBLE,0,52,comm);
    }
    else {
      start_num=0;
      for(i=0;i<n_proc;++i) {
        // c) receive
        if(i!=0) {
          MPI_Recv(buffer_rec_double,number_points_process[i],
                   MPI_DOUBLE,i,52,comm,&status);       
        }
        else {
          for(num=0;num<number_points_process[i];++num) {
            buffer_rec_double[num] = buffer_send_double[num];
          }
        }
        // d) print in file
        if(i<n_proc) {
          for(num=0;num<number_points_process[i];++num) {
            *Datei << " " << buffer_rec_double[num]
                   << "\n";
          }
        }
        start_num = start_num + number_points_process[i];
      }
    }


    // Teil 7: end remarks
    if(my_rank==0) {
      *Datei << "object " << '"' 
             << "irregular positions irregular connections"
             << '"' << " class field" << 
        endl;
      *Datei << "component " << '"' << "positions" << '"' << " value 1" << endl;
      *Datei << "component " 
             << '"' << "connections" << '"' << " value 2" << endl;
      *Datei << "component " << '"' << "data" << '"' << " value 3" << endl;
      *Datei << "end" << endl; 
    }      

    // Teil 8:
    if(I_am_active()) {
       delete(number_cells_process);
       delete(number_points_process);
    }
    delete(buffer_send_int);
    delete(buffer_send_double);
    if(my_rank==0) {
      delete(buffer_rec_int);
      delete(buffer_rec_double);
    }

    
  }



  MPI_Barrier(comm);
}

---