performance Optimierung

program/main.cc

@@ -98,12 +98,25 @@ int main(int argc, char** argv)
   u      = Cos(X*M_PI)*Cos(Y*M_PI)*Cos(Z*M_PI);
   
   
+  /*
   Interpolate_from_cell_to_point<Stencil_constant_z> interpolGrid(iteration);
+  interpolGrid.init(&blockgrid); 
+  interpolGrid.doInterpolation(&u_cell,&interpol);
+*/
 
-  interpolGrid.init(&blockgrid);
+  FastInterpolate_from_cell_to_point interpolGrid;
+  interpolGrid.init(&blockgrid); 
+  interpolGrid.doInterpolation(&u_cell,&interpol);
   
- interpolGrid.doInterpolation(&u_cell,&interpol);
+  
+    std::ofstream DATEIA;
 
+    DATEIA.open("interpol.vtk");
+    interpol.Print_VTK(DATEIA);
+    DATEIA.close();     
+     
+//     cout << " Maximum: " << Maximum(interpol) << " Minimum: " << Minimum(interpol) << endl;  
+  
   cout << " Maximum : " << L_infty(u) << endl; 
   cout << " Error : " << L_infty(interpol-u) << endl;
   

program/program.pro

@@ -17,7 +17,7 @@ unix: MOC_DIR = MOC_FILES
 
 
 # getting pathes:
-include(../config.pri)
+include(../configMy.pri)
 
 INCLUDEPATH += $${UGBLOCKS_PATH}
 

program/source/extemp/variable.cc

@@ -60,6 +60,24 @@ void Variable< std::complex<double> >::Update_back(int id_hex) const {
                               my_rank, comm);
 }
 
+template <>
+void Variable<double>::AddUpFromHex(int id_hex) const {
+    AddUp_back_fromHex<double>(id_hex, blockgrid,
+                              data_hexahedra, data_quadrangles,
+                              data_edges, data_points,
+                              my_rank, comm);
+}
+
+template <>
+void Variable< std::complex<double> >::AddUpFromHex(int id_hex) const {
+    AddUp_back_fromHex< std::complex<double> >(id_hex, blockgrid,
+                              data_hexahedra, data_quadrangles,
+                              data_edges, data_points,
+                              my_rank, comm);
+}
+
+
+
 template <>
 template <>
 void Variable<double>::Update<hexahedronEl>(int id_hex) const{

program/source/extemp/variable.h

@@ -306,13 +306,21 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
 		inline DTyp Give_data ( params_in, D3vector &posVector ) const ;
 		
 		void UpdateHexahedra(); ///> Updated alle Randdaten in allen Hexahedra
-        void UpdateHexahedraBack();
+                void UpdateHexahedraBack();
 		/// todo update nur falls notwendig
 		template <elementTyp TYP_EL>
 		void Update ( int id ) const;
 		
 		void Update_back(int id) const; ///> fuer interpolation slice 2D -> 3D
 				
+				
+                void AddUpFromHex();             ///> fuer Interpolation cell -> point
+                void AddUpFromHex(int id) const; ///> fuer Interpolation cell -> point	
+                void SetAllZero();               ///> fuer Interpolation cell -> point	                
+                void AddFromOne();               ///> fuer Interpolation cell -> point	
+	        template <class Zell>
+	        void AddFromZell(Zell& ZellVar); ///> fuer Interpolation cell -> point	
+				
 		
 		inline DTyp Give_data_quadrangle_middle ( int id, int i, int j, int Nx ) const;
 		inline DTyp Give_data_quadrangle_interior ( int id, int i, int j, int Nx ) const;
@@ -430,6 +438,15 @@ void Variable<DTyp>::UpdateHexahedraBack() {
      }
 }
 
+template <class DTyp>
+void Variable<DTyp>::AddUpFromHex() {
+     for(int id_hex=0;id_hex<ug->Give_number_hexahedra();++id_hex) {
+         AddUpFromHex(id_hex);
+     }
+}
+
+
+
 template <class DTyp>
 int Variable<DTyp>::checkThreshold ( double threshold, int kz ) {
 			int id_hex = 1;
@@ -1494,11 +1511,123 @@ void Variable<DTyp>::operator= ( const Expr<A>& a ) {
 }
 
 
+template <class DTyp>
+void Variable<DTyp>::SetAllZero() {
+	int Nx, Ny, Nz;
+	
+	bool useOpenMP_here = UGBlocks::useOpenMP;	
+	
+	// hexahedra	
+	for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {		
+			Nx = blockgrid->Give_Nx_hexahedron ( id );
+			Ny = blockgrid->Give_Ny_hexahedron ( id );
+			Nz = blockgrid->Give_Nz_hexahedron ( id );
+
+                        #pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(useOpenMP_here)
+			for(int k = 0;k <= Nz;++k ) {
+				for(int j = 0;j <= Ny;++j )
+				    for(int i = 0;i <= Nx;++i )
+					data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] = 0.0;
+			}
+		}
+	}
+	
+	// quadrangles
+	for ( int id = 0;id < ug->Give_number_quadrangles();++id ) {
+			
+		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {		
+			Nx = blockgrid->Give_Nx_quadrangle ( id );
+			Ny = blockgrid->Give_Ny_quadrangle ( id );
+
+                        // #pragma omp parallel for private(i) num_threads(UGBlocks::numThreadsToTake) if(useOpenMP_here)			
+			for(int j = 0;j <= Ny;++j )
+				for (int i = 0;i <= Nx;++i ) {
+					data_quadrangles[id][i+ ( Nx+1 ) *j] = 0.0;
+				}
+		}
+	}
+	
+	// edges
+	for ( int id = 0;id < ug->Give_number_edges();++id ) {
+		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
+			Nx = blockgrid->Give_Nx_edge ( id );
+			
+			for(int i = 0;i <= Nx;++i )
+				data_edges[id][i] = 0.0;
+		}
+	}
+	
+	// points
+	for(int id = 0;id < ug->Give_number_points();++id ) {
+		if(ug->Give_point ( id )->my_object ( my_rank ) ) {
+			data_points[id][0] = 0.0;
+		}
+	}
+}
 
 
 
 
+template <class DTyp>
+template <class Zell>
+void Variable<DTyp>::AddFromZell(Zell& ZellVar) {
+	int Nx, Ny, Nz;
+		
+	// hexahedra	
+	for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+	    if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {		
+			Nx = blockgrid->Give_Nx_hexahedron ( id );
+			Ny = blockgrid->Give_Ny_hexahedron ( id );
+			Nz = blockgrid->Give_Nz_hexahedron ( id );
+			
+			for(int k = 0;k < Nz;++k ) {
+				for(int j = 0;j < Ny;++j )
+				    for(int i = 0;i < Nx;++i ) {
+				        DTyp value = ZellVar.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
+					data_hexahedra[id][i   + (Nx+1) * (j   + (Ny+1)*k     ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j   + (Ny+1)*k     ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j+1 + (Ny+1)*k     ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*k     ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j   + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j   + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+				    }
+			}
+	    }
+	}
+}
 
+template <class DTyp>
+void Variable<DTyp>::AddFromOne() {
+	int Nx, Ny, Nz;
+		
+	// hexahedra	
+	for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+	    if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {		
+			Nx = blockgrid->Give_Nx_hexahedron ( id );
+			Ny = blockgrid->Give_Ny_hexahedron ( id );
+			Nz = blockgrid->Give_Nz_hexahedron ( id );
+			
+			
+			for(int k = 0;k < Nz;++k ) {
+				for(int j = 0;j < Ny;++j )
+				    for(int i = 0;i < Nx;++i ) {
+				        DTyp value = 1.0;
+					data_hexahedra[id][i   + (Nx+1) * (j   + (Ny+1)*k   ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j   + (Ny+1)*k   ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j+1 + (Ny+1)*k   ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*k   ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j   + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j   + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i   + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+					data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+				    }
+			}
+	    }
+	}
+}
 
 /////////////////////////////////////////////////////////////////////////
 // 5. Nconst , ... funtions for copy of data

program/source/extemp/variable_cc.h

@@ -1083,6 +1083,167 @@ assert(	rank_rec == 0);
 
 
 
+template <class DTyp>
+void AddUp_back_fromHex(int id_hex, Blockgrid* blockgrid, DTyp** data_hexahedra,
+                            DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
+                            int my_rank, MPI_Comm comm ) {
+	int i,j;
+	int id_quad, id_edge, id_point;
+	int Ni,Nj,Nconst;
+	int Nxedge;
+	int Nxquad,Nyquad;
+	int Nxhex,Nyhex,Nzhex;
+	int id_SW, id_SE, id_NW, idsp;
+	int id_W, id_E;
+	dir3D_sons c_SW ( ( dir3D_sons ) 0 ), c_SE ( ( dir3D_sons ) 0 ), c_NW ( ( dir3D_sons ) 0 );
+	dir3D_sons c_W ( ( dir3D_sons ) 0 ), c_E ( ( dir3D_sons ) 0 );
+	dir3D   c_x, c_y;
+	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();
+
+	int rank_send, rank_rec;  // for parallel
+
+	Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
+	Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
+	Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );
+
+	rank_rec = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel
+
+// keine MPI Paralleliserung!!!!!!!!!!!!!!
+assert(	rank_rec == 0);
+	
+	// quadrangles
+	for(int fc=0; fc<6; ++fc ) {
+	    id_quad = ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc );
+	    rank_send = ug->Give_quadrangle ( id_quad )->Give_rank(); // for parallel
+	    if(rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+	       Nxquad = blockgrid->Give_Nx_quadrangle ( id_quad );
+	       Nyquad = blockgrid->Give_Ny_quadrangle ( id_quad );
+
+	       // a) coordinate points of quadrangle
+	       id_SW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
+	       id_SE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
+	       id_NW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );
+
+	       // b) Find corresponding points in hexahedron
+	       for(int c=0; c<8; ++c ) {
+		   idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+		   if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
+		   else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
+		   else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
+	       }
+
+	       c_x = direction_from_to ( c_SW,c_SE );
+	       c_y = direction_from_to ( c_SW,c_NW );
+
+	       // x - direction
+	       Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+	       Ni     = Give_Ni (    ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+
+	       // y - direction
+	       Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
+	       Nj     = Give_Ni (    ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );
+
+	       // empty - direction
+	       Nconst = Give_Nconst ( opposite3D ( ( dir3D ) fc ),
+			                       Nxhex, Nyhex, Nzhex ) + Nconst;
+	       if(rank_send == rank_rec ) {    // for parallel
+		  for(j=1; j<Nyquad; ++j )
+		      for(i=1;i<Nxquad; ++i ) {
+			  // vorwaerts:
+			  //  data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
+			  //			   data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j];
+			  data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] +=
+                                          data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
+		      }
+	       }
+	       else {
+		  assert(false); 
+	       }
+	    }
+	}
+
+	// edges
+	for ( int ed=0; ed<12; ++ed ) {
+		id_edge = ug->Give_hexahedron ( id_hex )->Give_id_edge ( ( Edges_cell ) ed );
+		rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
+		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+			Nxedge = blockgrid->Give_Nx_edge ( id_edge );
+
+			// a) coordinate points of edge
+			id_W = ug->Give_edge ( id_edge )->Give_id_corner_W();
+			id_E = ug->Give_edge ( id_edge )->Give_id_corner_E();
+
+			// b) Find corresponding points in hexahedron
+			for ( int c=0; c<8; ++c ) {
+				idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+				if ( idsp==id_W ) c_W = ( dir3D_sons ) c;
+				else if ( idsp==id_E ) c_E = ( dir3D_sons ) c;
+			}
+
+			c_x = direction_from_to ( c_W,c_E );
+
+			// x - direction
+			Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+			Ni     = Give_Ni    ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+
+			// A. empty - direction
+			Nconst = Give_Nconst ( opposite3D ( OrthoA ( ( Edges_cell ) ed ) ),
+			                       Nxhex, Nyhex, Nzhex ) + Nconst;
+
+			// B. empty - direction
+			Nconst = Give_Nconst ( opposite3D ( OrthoB ( ( Edges_cell ) ed ) ),
+			                       Nxhex, Nyhex, Nzhex ) + Nconst;
+
+			if(rank_send == rank_rec ) {    // for parallel
+			   for(i=1; i<Nxedge; ++i ) {
+			       // vorwaerts
+			       //data_hexahedra[id_hex][i*Ni+Nconst] = data_edges[id_edge][i];
+                               data_edges[id_edge][i] += data_hexahedra[id_hex][i*Ni+Nconst];
+			   }
+			} 
+			else {
+			   assert(false); 
+			}
+		}
+	}
+
+	// points
+	for ( int c=0; c<8; ++c ) {
+		id_point = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+		rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
+		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+
+			// x- direction
+			Nconst = Give_Nconst ( opposite3D ( xCoord ( ( dir3D_sons ) c ) ),
+			                       Nxhex, Nyhex, Nzhex );
+
+			// y- direction
+			Nconst = Give_Nconst ( opposite3D ( yCoord ( ( dir3D_sons ) c ) ),
+			                       Nxhex, Nyhex, Nzhex ) + Nconst;
+
+			// z- direction
+			Nconst = Give_Nconst ( opposite3D ( zCoord ( ( dir3D_sons ) c ) ),
+			                       Nxhex, Nyhex, Nzhex ) + Nconst;
+
+			if(rank_send == rank_rec )      // for parallel
+                           // vorwaerts
+			   //data_hexahedra[id_hex][Nconst] = data_points[id_point][0];
+              {
+                data_points[id_point][0] += data_hexahedra[id_hex][Nconst];
+            }
+			else {
+			  assert(false);
+			}
+		}
+	}
+}
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+
 
 
 

program/source/interpol/interpolCellVar.cc

@@ -26,6 +26,50 @@ using namespace ::_COLSAMM_;
 #include "assert.h"
 
 
+FastInterpolate_from_cell_to_point::FastInterpolate_from_cell_to_point() {
+     weight = NULL;
+}
+
+void FastInterpolate_from_cell_to_point::deleteData() {
+   if(weight!=NULL) delete weight;
+   weight = NULL;
+}
+
+void FastInterpolate_from_cell_to_point::init(Blockgrid* blockgrid) {
+     deleteData();
+     weight = new Variable<double>(*blockgrid);
+     weight->SetAllZero();
+     weight->AddFromOne();
+     weight->AddUpFromHex();
+     weight->UpdateHexahedra();
+     
+/*     
+    std::ofstream DATEIA;
+
+    DATEIA.open("weight.vtk");
+    weight->Print_VTK(DATEIA);
+    DATEIA.close();     
+     
+     cout << " Maximum: " << Maximum(*weight) << " Minimum: " << Minimum(*weight) << endl;
+     */
+}
+
+void FastInterpolate_from_cell_to_point::doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint) {
+     assert(weight!=NULL);
+     varPoint->SetAllZero();
+     varPoint->AddFromZell(*varCell);
+     varPoint->AddUpFromHex();
+     
+     (*varPoint) = (*varPoint) / (*weight);
+     varPoint->UpdateHexahedra();     
+     
+     
+     
+}
+
+
+/////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////
 
 
 template<>
@@ -68,3 +112,4 @@ void Interpolate_from_cell_to_point<Stencil_constant_z>::init(Blockgrid* blockgr
 	
      }
 }
+

program/source/interpol/interpolCellVar.h

@@ -48,6 +48,27 @@
 #define INTERPOLCELLVAR_H
 
 
+/** \addtogroup InterpolationOperators **/
+/* @{ */  
+/**
+ * Interpolation from Cell_variable to Variable <br>
+ * FAST VERSION <br>
+ * Sometimes this version is a little bit less accurate than Interpolate_from_cell_to_point in case of smooth functions <br>
+ * But in case of a non smooth function, this version is more appropriate
+**/
+class FastInterpolate_from_cell_to_point {
+ public:
+  FastInterpolate_from_cell_to_point();
+
+  void init(Blockgrid* blockgrid_);    
+  void deleteData();
+  void doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint);
+  
+  ~FastInterpolate_from_cell_to_point() { deleteData(); }
+private:
+  Variable<double>* weight;  
+};
+
 
 
 
@@ -56,8 +77,11 @@
 
 
 /**
- * interpolation by Gauss-Seidel from Cell_variable to Variable
- * template might be Stencil_constant_z
+ * Interpolation by Gauss-Seidel from Cell_variable to Variable <br>
+ * template might be Stencil_constant_z <br>
+ * SLOW VERSION <br>
+ * Sometimes this version is more accurate than FastInterpolate_from_cell_to_point <br>
+ * An example showed an increase of factor 3 of accuracy in case of a smooth function 
 **/
 template<class StenTyp>
 class Interpolate_from_cell_to_point {
@@ -81,6 +105,8 @@ private:
   
   int numberGSinterations;
 };
+/* @} */
+
 
 template<class StenTyp>
 void Interpolate_from_cell_to_point<StenTyp>::doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint) {
@@ -125,4 +151,5 @@ void Interpolate_from_cell_to_point<StenTyp>::deleteData() {
      }
 }
 
+
 #endif // INTERPOLCELLVAR_H