/*============================================================================
 * This function is called at the end of each time step, and has a very
 *  general purpose
 *  (i.e. anything that does not have another dedicated user function)
 *============================================================================*/

/* Code_Saturne version 5.0.5 */

/*
  This file is part of Code_Saturne, a general-purpose CFD tool.

  Copyright (C) 1998-2017 EDF S.A.

  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.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

/*----------------------------------------------------------------------------*/

#include "cs_defs.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>
#include <stdio.h>

#if defined(HAVE_MPI)
#include <mpi.h>
#endif

/*----------------------------------------------------------------------------
 * PLE library headers
 *----------------------------------------------------------------------------*/

#include <ple_coupling.h>

/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/

#include "bft_mem.h"
#include "bft_error.h"
#include "bft_printf.h"

#include "cs_base.h"
#include "cs_field.h"
#include "cs_field_pointer.h"
#include "cs_field_operator.h"
#include "cs_math.h"
#include "cs_mesh.h"
#include "cs_mesh_quantities.h"
#include "cs_halo.h"
#include "cs_halo_perio.h"
#include "cs_log.h"
#include "cs_parall.h"
#include "cs_parameters.h"
#include "cs_physical_constants.h"
#include "cs_post_util.h"
#include "cs_rotation.h"
#include "cs_time_moment.h"
#include "cs_time_step.h"
#include "cs_turbomachinery.h"
#include "cs_selector.h"

#include "cs_post.h"

/*----------------------------------------------------------------------------
 *  Header for the current file
 *----------------------------------------------------------------------------*/

#include "cs_prototypes.h"
#include "cs_base.h"
#include "cs_mesh.h"
#include "cs_mesh_quantities.h"
#include "cs_mesh_bad_cells.h"
#include "cs_probe.h"
#include "cs_volume_zone.h"
#include "cs_domain.h"
/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS

/*----------------------------------------------------------------------------*/
/*!
 * \file cs_user_extra_operations-turbomachinery.c
 *
 * \brief This function is called at the end of each time step, and has a very
 * general purpose (i.e. anything that does not have another dedicated
 * user function).
 *
 * This example is a part of the \subpage turbomachinery example.
 */
/*----------------------------------------------------------------------------*/

/*============================================================================
 * User function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------
 * User function to express a vector in the cyclindrical coordinates system of
 * a given rotation zone.
 *
 * parameters:
 *   r      <-- pointer to rotation structure
 *   coords <-- coordinates at point
 *   v      <-- vector components in cartesian coordinates system
 *   vc     --> vector components in cylindrical coordinates system
 *----------------------------------------------------------------------------*/
 

/*----------------------------------------------------------------------------
 * User function to find the closest cell and the associated rank of a point
 * in a given rotation zone.
 *
 * parameters:
 *   r       <-- id of the rotation zone
 *   coords  <-- point coordinates
 *   node    --> cell id
 *   rank    --> rank id
 *----------------------------------------------------------------------------*/


static void
_findpt_r(const cs_real_3_t     coords,
          cs_lnum_t            *node,
          cs_lnum_t            *rank)
{
  cs_real_t d[3];
  cs_lnum_t n_cells = cs_glob_mesh-> n_cells;
  cs_real_3_t *cell_cen = (cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;


  *node = (int)(1);
  for (int i = 0; i < 3; i++){
    d[i] = coords[i] - cell_cen[*node][i];
  }
  cs_real_t dis2mn = cs_math_3_square_norm(d);

  for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
      for (int i = 0; i < 3; i++){
        d[i] = coords[i] - cell_cen[cell_id][i];
	  }
	  cs_real_t dis2 = cs_math_3_square_norm(d);
	  if (dis2 < dis2mn) {
		*node = cell_id;
		dis2mn = dis2;
		}
	}
  cs_parall_min_id_rank_r(node, rank, dis2mn);
}

void
cs_user_extra_operations(void)
{

  /* Mesh-related variables */

  const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
  double d;
  const cs_real_3_t  *restrict cell_cen
    = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->cell_cen;

  /* 0. Initialization
     ================= */

  /* Local constants defintion */

  const cs_real_t _PI = acos(-1.);
  const cs_real_t _G = 9.81;

  /* Example 2: extraction of a velocity profile in cylindrical corordinates
     ======================================================================= */

  if (cs_glob_time_step->nt_cur == cs_glob_time_step->nt_max){
	cs_lnum_t n_cells = cs_glob_mesh-> n_cells;
    cs_real_3_t *vel = (cs_real_3_t *)CS_F_(u)->val;
	double l;
	double u, v, w;
	double mean_u, mean_v, mean_w;
    FILE *f1 = NULL, *f2 = NULL, *f3 = NULL;
	int j=0;
	double sum_u, sum_v, sum_w;
    /* Only process of rank 0 (parallel) or -1 (scalar) writes to this file. */
	d = (double)n_cells;
    f1 = fopen("sum.dat","w");

    f2 = fopen("mean.dat","w");
	
	f3 = fopen("z_dir.dat","w");

    cs_real_t radius = 0.00805;
    cs_real_t axicoo = 0.585;
    cs_lnum_t npoint = 100;

	for (double radius_curr = 0.0060; radius_curr<radius; radius_curr+=0.0001){
		j = 0;
		sum_u = 0.0;
		sum_v = 0.0;
		sum_w = 0.0;
		for (double z_dir = 0.000; z_dir<axicoo; z_dir+=0.00585){
			for (cs_lnum_t point_id = 0; point_id < npoint; point_id++) {
			  cs_real_t xth0 = -_PI/npoint;
			  cs_real_3_t xyz = {radius_curr*cos((double)point_id/npoint*2.*_PI),
								 radius_curr*sin((double)point_id/npoint*2.*_PI),
								 z_dir};
			  cs_lnum_t cell_id, rank_id;
			  int a1 = (int)rank_id;
			  _findpt_r(xyz, &cell_id, &rank_id);
			  /*findpt(n_cells, n_cells, cell_cen, xyz[0], xyz[1], xyz[2], &cell_id, &rank_id);*/
			  if(a1==(int)rank_id){
				  j++;
				  u = vel[cell_id][0];
				  v = vel[cell_id][1];
				  w = vel[cell_id][2];
				  sum_u += (double)u;
				  sum_v += (double)v;
			      sum_w += (double)w;
				  l = sqrt((double)cell_cen[cell_id][0]*(double)cell_cen[cell_id][0] + (double)cell_cen[cell_id][1]*(double)cell_cen[cell_id][1]);
				  fprintf(f3," %e\t%e \t %e\t%e\t %e\t%d \t %e\n", cell_cen[cell_id][0], cell_cen[cell_id][1], cell_cen[cell_id][2],  z_dir, radius_curr, (int)rank_id, l);
				  if(radius_curr-0.0003<l && l<radius_curr+0.0003 && cell_cen[cell_id][2] < z_dir+0.01 && cell_cen[cell_id][2] >z_dir-0.01){
				}
			  }
			}	/* end of loop on point_id */
		}
		mean_u = (double)sum_u /j;
		mean_v = (double)sum_v /j;
		mean_w = (double)sum_w /j;
		fprintf(f1,"  %17.9e%17.9e%17.9e%17.9e\t%d\n", radius_curr, sum_u, sum_v, sum_w, j);
		fprintf(f2,"  %17.9e%17.9e%17.9e\n", mean_u, mean_v, mean_w);
	}
    fclose(f1);
    fclose(f2);
	fclose(f3);
  } /* end of if statement on current time step */
}

END_C_DECLS