#include "precice/SolverInterfaceC.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <petsc.h>

typedef struct {
	PetscScalar x, y, z;
} Cmpnts;

typedef struct {
	int		nbnumber;
	int		n_v, n_elmt;	// number of vertices and number of elements
	int		my_n_v, my_n_elmt;	// seokkoo, my proc
	int		*nv1, *nv2, *nv3;	// Node index of each triangle
	PetscReal	*nf_x, *nf_y, *nf_z;	// Normal direction
	//PetscReal	*nf_x0, *nf_y0, *nf_z0;	// Normal direction
	PetscReal	*x_bp, *y_bp, *z_bp;	// Coordinates of IBM surface nodes
	PetscReal	*x_bp0, *y_bp0, *z_bp0;
    PetscReal       *x_bp_i, *y_bp_i, *z_bp_i;
	PetscReal	*x_bp_o, *y_bp_o, *z_bp_o;

	//PetscReal	x_bp_in[101][3270], y_bp_in[101][3270], z_bp_in[101][3270];
	Cmpnts		*u, *uold, *urm1, *u_pre;
   	// Acceleration


	PetscReal     	*x_bp_pre, *y_bp_pre, *z_bp_pre;
	PetscReal     	*error;
	PetscReal     	*error_force;

	// Aitken
	PetscReal   	Aitken_factor;
	PetscReal   	Aitken_factor_force;

	// shear stress projection
	PetscReal       *Txx,*Txy, *Txz, *Tyy, *Tzz, *Tyz;
  	PetscReal       *Traction_x, *Traction_y, *Traction_z;
  	PetscReal       *Traction_x_pre, *Traction_y_pre, *Traction_z_pre;
  	PetscReal       *Area_x, *Area_y, *Area_z;
  	PetscReal      	*SumP_Tx,*SumP_Ty,*SumP_Tz,*AbsP_tau;
 	PetscReal     	*SumM_Tx,*SumM_Ty,*SumM_Tz,*AbsM_tau;

	// Added 4/1/06 iman
	PetscReal     	*dA ;         // area of an element
	// Added 4/2/06 iman
	PetscReal     	*nt_x, *nt_y, *nt_z; //tangent dir
	PetscReal     	*ns_x, *ns_y, *ns_z; //azimuthal dir
	// Added 6/4/06 iman
	//Cmpnts        *cent; // center of the element
	PetscReal     	*cent_x,*cent_y,*cent_z;

	// element center displacement
	PetscReal	*disp_x, *disp_y, *disp_z, *disp_theta;
	PetscReal	*ux_struct, *uy_struct, *uz_struct, *utheta_struct;

	// for radius check
	Cmpnts 		*qvec;
	PetscReal 	*radvec;

	//seokkoo

	//int *_nv1, *_nv2, *_nv3;	// for rank 0
	PetscReal 	*count;
	//int *local2global_elmt;
	//int total_n_elmt, total_n_v;

  //PetscReal *_x_bp, *_y_bp, *_z_bp;	// for rank 0
	PetscReal 	*shear;
	Cmpnts 		*rel_velocity;	// flow velocity - body velocity


        // add begin (xiaolei)
	// fsi_fem
	int 		bctype, thin, masscorrection;
  	int       	fsi_type; // 0 - stationary ; 1 -- prescribed motion ; 2 -- Solid body  ; 3-- Deformable body

	PetscReal 	*Tmprt_lagr, *Ftmprt_lagr, *tmprt;

	    PetscReal 	*F_lagr_x, *F_lagr_y, *F_lagr_z; // force at the IB surface points (lagrange points)
        PetscReal 	*dF_lagr_x, *dF_lagr_y, *dF_lagr_z;

        PetscReal 	*Ft_lagr_avg, *Fa_lagr_avg, *Fr_lagr_avg; // force at the IB surface points (lagrange points)
        PetscReal 	*U_lagr_x, *U_lagr_y, *U_lagr_z; // estimated fluid velocity at Lagrange points


	PetscReal 	*Urelmag; // relative incoming velocity for actuator model

	Cmpnts		*Urel; // vector of the relative incoming velocity
	Cmpnts		*Uinduced; // vector of the induced velocity
	Cmpnts		*circulation; // circulation vector on the blade
	Cmpnts		*liftdirection; // direction of the lift
	Cmpnts		*rotationdirection; // direction of the lift



	Cmpnts		*Urel_mean; // vector of the relative incoming velocity
	Cmpnts		*Uinduced_mean; // vector of the induced velocity
	Cmpnts		*circulation_mean; // circulation vector on the blade
	Cmpnts		*liftdirection_mean; // direction of the lift

        int        	*i_min, *i_max, *j_min, *j_max, *k_min, *k_max;

        /* ACL */
        PetscReal       *angle_attack, *angle_twist, *chord_blade; // twist angle and chord length at each grid point
        PetscReal       *CD, *CL;
	PetscReal	*pitch;
	PetscReal	*pitch_old;
	PetscReal	*pitch_IPC;
	PetscReal	pitch_min;

	PetscReal 	U_ref;



	PetscReal	*dhx, *dhy, *dhz;

	PetscReal 	CD_bluff;
	PetscReal 	friction_factor, pressure_factor;

	PetscReal	*frictionfactor;

	PetscReal 	axialforcecoefficient, tangentialforcecoefficient;
	PetscReal 	axialprojectedarea, tangentialprojectedarea;

	PetscReal	dh;

	PetscReal	indf_axis, Tipspeedratio, CT, indf_tangent;

	int 		num_CT; // number of input CT nodes
	PetscReal	rCTInp[200], CTInp[200]; // Thrust coeffcient from a file HARD CODING
	PetscReal	*CTdisk;


        PetscReal 	*Fr_mean, *Fa_mean, *Ft_mean; // force at the IB surface points (lagrange points)
        PetscReal 	*Ur_mean, *Ua_mean, *Ut_mean;
	PetscReal 	*AOA_mean, *Urelmag_mean;


        PetscReal 	*AOAAOA_mean, *FFa_mean, *FFt_mean;


	PetscReal	*centIP_x, *centIP_y, *centIP_z;

        int        	*iIP_min, *iIP_max, *jIP_min, *jIP_max, *kIP_min, *kIP_max;
        PetscReal 	*U_IPlagr_x, *U_IPlagr_y, *U_IPlagr_z, *P_IPlagr;

	PetscReal	*dh_IP;

	PetscReal	*Nut_lagr, *Shear_lagr_x, *Shear_lagr_y, *Shear_lagr_z;
	PetscReal	*ShearDesired_lagr_x, *ShearDesired_lagr_y, *ShearDesired_lagr_z;
	PetscReal	*UU_lagr_x, *UU_lagr_y, *UU_lagr_z;

	PetscReal	*random_color;

	int		num_cf;
	PetscReal	r_in[200], cf_in[200]; // Readed friction coefficient from a file  HARD CODING

	PetscReal	forcecoeffx_actuator, forcecoeffy_actuator, forcecoeffz_actuator, Uref_actuator, Lref_actuator;


	int		num_circulationInp;
	PetscReal	*circulationInp, *r_circulationInp; // simulation with circulation specified along the blade
	PetscReal	*circulationSpecified;

	int		num_blade, num_foiltype;



	PetscReal	r_nearhubinflowcorrection;

	PetscReal	coeff_SG;

	PetscReal	r_nacelle;

        // add end (xiaolei)

	// sediment begin
    	// Added 1/10/2010 ali
    	PetscReal     	*deltz_p_us, *deltz_p_ds, *A_us, *A_ds;
    	PetscReal     	*BShS, *BCShS, *Shvel;
    	PetscReal     	*SCont;
    	PetscReal     	*C;
    	PetscReal     	*max_bed_angle, *SedFlux;
    	int      	*Rigidity, *Mobility;
    	PetscReal     	*Dely, *cent_y_old, *y_bp_l;
    	PetscReal     	*Delz, *cent_z_old, *z_bp_l;
    	PetscReal      	Deltab;
    	PetscReal     	*elmt_depth;
    	PetscReal     	*netflux_old;
    	PetscReal     	*F_flux_12,*F_flux_13,*F_flux_23;
    	PetscReal     	dtime,time_bedchange;
    	PetscReal 	*cent_zl, *cent_z_AVE;



	int 		*node_boundary;
	int 		*elmt_boundary;

	// xiaolei add SEDI

    	int  		*color;
    	int 		*s2l; // actuator line element index for each actuator surface element

	// end add SEDI

	// sediment end
	//
	// fsi_IBM
	 PetscReal      F_x, F_y, F_z;

  	// Total surface area
     	PetscReal      	Area;

} IBMNodes;

typedef struct {
	PetscReal    x_c,y_c,z_c; // center of rotation(mass)
    PetscReal    x_c0,y_c0,z_c0;
}FSInfo;

typedef struct {
	PetscInt dim;
  PetscInt writeDataID;
  PetscInt readDataID;
  PetscInt meshID;
  PetscInt numberOfVertices;
  PetscReal *vertices;
  PetscReal *readData;
  PetscReal *writeData;
  PetscInt *vertexIDs;
  const char *meshName;
  const char *writeDataName;
  const char *readDataName;
  const char *location;

} PreciceInfo;

//wzw: add for output the surface of the airfoil:velocity
PetscErrorCode velocity_surface_out(IBMNodes *ibm, int ti, int ibi)
{
  int rank,i;
  int time_intv=100;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  if (!rank) {
    if ((ti%time_intv) == 0) {
    FILE *f;
    char filen[80];
    sprintf(filen, "VelocityInterface%06d_%2.2d.dat",ti,ibi+1);
    f = fopen(filen, "w");
	    
    int N_block = 100 * 1024 * 1024; // 100Mb
    char str[256];
    char carriage_return = '\n';
    //std::vector<char> large_buffer(N_block); // Initialize the large_buffer vector with the specified block size
	    
      PetscFPrintf(PETSC_COMM_WORLD, f, "Variables=x,y,z,u_x,u_y,u_z,n_x,n_y,n_z,nt_x,nt_y,nt_z\n");
      PetscFPrintf(PETSC_COMM_WORLD, f, "ZONE T='TRIANGLES', N=%d, E=%d, F=FEBLOCK, ET=TRIANGLE, VARLOCATION=([1-3]=NODAL,[4-12]=CELLCENTERED)\n", ibm->n_v, ibm->n_elmt);

      for (i=0; i<ibm->n_v; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->x_bp[i]);
      }
      for (i=0; i<ibm->n_v; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->y_bp[i]);
      }
      for (i=0; i<ibm->n_v; i++) {	
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->z_bp[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->U_lagr_x[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->U_lagr_y[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->U_lagr_z[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_x[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_y[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_z[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_x[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_y[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_z[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%d %d %d\n", ibm->nv1[i]+1, ibm->nv2[i]+1, ibm->nv3[i]+1);
      }
      fclose(f);
    } 
  }
  	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: velocity_surface_out  \n");

  return(0);
}

PetscErrorCode Force_surface_out(IBMNodes *ibm, int ti, int ibi)
{
  int rank,i;
  int time_intv=100;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  
  char path[256];
  sprintf(path, ".");

  if (!rank) {
    if ((ti%time_intv) == 0) {
    FILE *f;
    char filen[80];
    sprintf(filen, "ForcesInterface%06d_%2.2d.dat",ti,ibi+1);
    f = fopen(filen, "w");
	    
    int N_block = 100 * 1024 * 1024; // 100Mb
    char str[256];
    char carriage_return = '\n';
    //std::vector<char> large_buffer(N_block); // Initialize the large_buffer vector with the specified block size
	    
      PetscFPrintf(PETSC_COMM_WORLD, f, "Variables=x,y,z,F_x,F_y,F_z,n_x,n_y,n_z,nt_x,nt_y,nt_z\n");
      PetscFPrintf(PETSC_COMM_WORLD, f, "ZONE T='TRIANGLES', N=%d, E=%d, F=FEBLOCK, ET=TRIANGLE, VARLOCATION=([1-3]=NODAL,[4-12]=CELLCENTERED)\n", ibm->n_v, ibm->n_elmt);

      for (i=0; i<ibm->n_v; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->x_bp[i]);
      }
      for (i=0; i<ibm->n_v; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->y_bp[i]);
      }
      for (i=0; i<ibm->n_v; i++) {	
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->z_bp[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->F_lagr_x[i]/ibm->dA[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->F_lagr_y[i]/ibm->dA[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->F_lagr_z[i]/ibm->dA[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_x[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_y[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nf_z[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_x[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_y[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%e\n", ibm->nt_z[i]);
      }
      for (i=0; i<ibm->n_elmt; i++) {
	PetscFPrintf(PETSC_COMM_WORLD, f, "%d %d %d\n", ibm->nv1[i]+1, ibm->nv2[i]+1, ibm->nv3[i]+1);
      }
      fclose(f);
    } 
  }
  	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: Force_surface_out  \n");

  return(0);
}

//wzw add for read surface data
PetscErrorCode Surface_read_ucd(IBMNodes *ibm, int ibi, FSInfo *fsi, char fname[80], double reflength)
{
  	int	rank;
  	int	n_v , n_elmt ;
  	PetscReal	*x_bp , *y_bp , *z_bp ;
  	int	*nv1 , *nv2 , *nv3 ;
  	PetscReal	*nf_x, *nf_y, *nf_z;
  	int	i,ii;
  	int	n1e, n2e, n3e;
  	PetscReal	dx12, dy12, dz12, dx13, dy13, dz13;
  	PetscReal     dr;
  	PetscReal     *dA ;//area
  	PetscReal	*nt_x, *nt_y, *nt_z;
  	PetscReal	*ns_x, *ns_y, *ns_z;
    char path[256];
    sprintf(path, ".");

  	char   ss[20];
  	char string[128];

      	//PetscPrintf(PETSC_COMM_WORLD, "xc=%le, yc=%le, zc=%le for %i th turbine\n", fsi->x_c, fsi->y_c, fsi->z_c, ibi);
      	//PetscPrintf(PETSC_COMM_WORLD, "nx=%le, ny=%le, nz=%le for %i th turbine\n", fsi->nx_tb, fsi->ny_tb, fsi->nz_tb, ibi);

  	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  	if(!rank) { // root processor read in the data
    		FILE *fd;
    		char filen[160];
    		sprintf(filen,"%s/%s" , path, fname);
    		PetscPrintf(PETSC_COMM_SELF, "READ %s\n", filen);

    		fd = fopen(filen, "r");
    		if (!fd) printf("Cannot open %s !!", filen),exit(0);
    		else printf("Opened %s !\n", filen);
    		n_v =0;

    		if (fd) {
      			fgets(string, 128, fd);
      			fgets(string, 128, fd);
      			fgets(string, 128, fd);

      			fscanf(fd, "%i %i %i %i %i",&n_v,&n_elmt,&ii,&ii,&ii);
      			PetscPrintf(PETSC_COMM_SELF, "number of nodes & elements %d %d\n",n_v, n_elmt);

      			ibm->n_v = n_v;
      			ibm->n_elmt = n_elmt;

      			MPI_Bcast(&(ibm->n_v), 1, MPI_INT, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(&(ibm->n_elmt), 1, MPI_INT, 0, PETSC_COMM_WORLD);

      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp));	// added by seokkoo 03.04.2009
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp));

      			x_bp = ibm->x_bp;	// seokkoo
      			y_bp = ibm->y_bp;
      			z_bp = ibm->z_bp;

      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_i));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_i));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_i));

      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_o));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_o));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_o));

      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp0));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp0));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp0));



      			PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->u));
      			PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->uold));
      			PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->urm1));

			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->tmprt));



		      	PetscMalloc(n_v*sizeof(int), &ibm->node_boundary); // xiaolei add
      			PetscMalloc(n_elmt*sizeof(int), &ibm->elmt_boundary); // xiaolei add

     	 		PetscMalloc(3*n_elmt*sizeof(PetscReal), &(ibm->error_force)); // xiaolei add


			// fsi_fem
      			PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->u_pre));

      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_pre));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_pre));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_pre));

			// Error Vector
			PetscMalloc(9*n_v*sizeof(PetscReal), &(ibm->error));

		      	PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_x);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_y);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_z);

      			// xiaolei add
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_x_pre);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_y_pre);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_z_pre);

	      		// xiaolei added
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Tx);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Ty);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Tz);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->AbsP_tau);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Tx);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Ty);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Tz);
			PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->AbsM_tau);



      			for (i=0; i<n_v; i++) {
				fscanf(fd, "%i %le %le %le", &ii, &x_bp[i], &y_bp[i], &z_bp[i]);//, &t, &t, &t);

				x_bp[i]=x_bp[i]/reflength;
				y_bp[i]=y_bp[i]/reflength;
				z_bp[i]=z_bp[i]/reflength;

			        ibm->x_bp_i[i] = x_bp[i];
			        ibm->y_bp_i[i] = y_bp[i];
			        ibm->z_bp_i[i] = z_bp[i];

			        x_bp[i] += fsi->x_c;
			        y_bp[i] += fsi->y_c;
			        z_bp[i] += fsi->z_c;

			        ibm->x_bp0[i] = x_bp[i];
			        ibm->y_bp0[i] = y_bp[i];
			        ibm->z_bp0[i] = z_bp[i];


				ibm->x_bp[i] = x_bp[i];
				ibm->y_bp[i] = y_bp[i];
				ibm->z_bp[i] = z_bp[i];

				ibm->x_bp_o[i] = x_bp[i];
				ibm->y_bp_o[i] = y_bp[i];
				ibm->z_bp_o[i] = z_bp[i];

				ibm->tmprt[i] = 0.;

				ibm->u[i].x = 0.;
				ibm->u[i].y = 0.;
				ibm->u[i].z = 0.;

				ibm->uold[i].x = 0.;
				ibm->uold[i].y = 0.;
				ibm->uold[i].z = 0.;

				ibm->urm1[i].x = 0.;
				ibm->urm1[i].y = 0.;
				ibm->urm1[i].z = 0.;
			}
      			i=0;
		      	PetscPrintf(PETSC_COMM_WORLD, "xyz_bp %le %le %le\n", x_bp[i], y_bp[i], z_bp[i]);

      			MPI_Bcast(ibm->x_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->y_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->z_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

      			MPI_Bcast(ibm->x_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->y_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->z_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

      			MPI_Bcast(ibm->x_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->y_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->z_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

		      	MPI_Bcast(ibm->x_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->y_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
      			MPI_Bcast(ibm->z_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);


      			PetscMalloc(n_elmt*sizeof(int), &nv1);
      			PetscMalloc(n_elmt*sizeof(int), &nv2);
      			PetscMalloc(n_elmt*sizeof(int), &nv3);

      			PetscMalloc(n_elmt*sizeof(PetscReal), &nf_x);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &nf_y);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &nf_z);

      			PetscMalloc(n_elmt*sizeof(PetscReal), &dA); //Area

      			PetscMalloc(n_elmt*sizeof(PetscReal), &nt_x);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &nt_y);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &nt_z);

      			PetscMalloc(n_elmt*sizeof(PetscReal), &ns_x);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ns_y);
      			PetscMalloc(n_elmt*sizeof(PetscReal), &ns_z);

      			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_x));
      			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_y));
      			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_z));

	      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_x));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_y));
      			PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_z));
      			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->disp_theta));


			// rotor model
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Urelmag_mean));

        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Fr_mean));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ft_mean));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Fa_mean));

	        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->AOAAOA_mean));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->FFt_mean));
	        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->FFa_mean));




        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ur_mean));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ut_mean));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ua_mean));

        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->AOA_mean));


        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Urelmag));

        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_x));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_y));
        		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_z));

	       		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_x));
	      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_y));
      			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_z));

		        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_x));
		        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_y));
		        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_z));


 		       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ftmprt_lagr));
			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Tmprt_lagr));


		        PetscMalloc(n_elmt*sizeof(int), &(ibm->i_min));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->j_min));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->k_min));

		        PetscMalloc(n_elmt*sizeof(int), &(ibm->i_max));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->j_max));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->k_max));

			PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhx));
       		        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhy));
                	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhz));


        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->angle_attack));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->angle_twist));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->chord_blade));

                       	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CD));
                       	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CL));


		        PetscMalloc(n_elmt*sizeof(int), &(ibm->iIP_min));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->jIP_min));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->kIP_min));

		        PetscMalloc(n_elmt*sizeof(int), &(ibm->iIP_max));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->jIP_max));
		        PetscMalloc(n_elmt*sizeof(int), &(ibm->kIP_max));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->random_color));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_x));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_y));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_z));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_x));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_y));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_z));

        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Nut_lagr));
        		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->dh_IP));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->P_IPlagr));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_x));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_y));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_z));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_x));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_y));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_z));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_x));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_y));
	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_z));

	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Urel));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Uinduced));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->circulation));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->liftdirection));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->rotationdirection));

	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Urel_mean));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Uinduced_mean));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->circulation_mean));
	        	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->liftdirection_mean));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->circulationSpecified));




	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->frictionfactor));

			PetscMalloc(n_elmt*sizeof(int), &(ibm->color));
			PetscMalloc(n_elmt*sizeof(int), &(ibm->s2l));

	        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CTdisk));

      			for (i=0; i<n_elmt; i++) {
				fscanf(fd, "%i %i %s %i %i %i\n", &ii,&ii, &ss, nv1+i, nv2+i, nv3+i);
				nv1[i] = nv1[i] - 1; nv2[i] = nv2[i]-1; nv3[i] = nv3[i] - 1;
		      	}
		        ibm->nv1 = nv1; ibm->nv2 = nv2; ibm->nv3 = nv3;

		        i=0;
		        PetscPrintf(PETSC_COMM_WORLD, "nv %d %d %d\n", nv1[i], nv2[i], nv3[i]);

		        fclose(fd);
 	      	}

	      	for (i=0; i<n_elmt; i++) {

			n1e = nv1[i]; n2e =nv2[i]; n3e = nv3[i];
		        dx12 = x_bp[n2e] - x_bp[n1e];
		        dy12 = y_bp[n2e] - y_bp[n1e];
		        dz12 = z_bp[n2e] - z_bp[n1e];

		        dx13 = x_bp[n3e] - x_bp[n1e];
		        dy13 = y_bp[n3e] - y_bp[n1e];
		        dz13 = z_bp[n3e] - z_bp[n1e];

		        nf_x[i] = dy12 * dz13 - dz12 * dy13;
		        nf_y[i] = -dx12 * dz13 + dz12 * dx13;
		        nf_z[i] = dx12 * dy13 - dy12 * dx13;

		        dr = sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i] + nf_z[i]*nf_z[i]);

		        nf_x[i] /=dr; nf_y[i]/=dr; nf_z[i]/=dr;

		      	if ((((1.-nf_z[i])<=1e-6 )&((-1.+nf_z[i])<1e-6))|
			  (((nf_z[i]+1.)<=1e-6 )&((-1.-nf_z[i])<1e-6))) {
				ns_x[i] = 1.;
				ns_y[i] = 0.;
				ns_z[i] = 0. ;

				nt_x[i] = 0.;
				nt_y[i] = 1.;
				nt_z[i] = 0.;
		        } else {
				ns_x[i] =  nf_y[i]/ sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i]);
				ns_y[i] = -nf_x[i]/ sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i]);
				ns_z[i] = 0. ;

				nt_x[i] = -nf_x[i]*nf_z[i]/ sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i]);
				nt_y[i] = -nf_y[i]*nf_z[i]/ sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i]);
				nt_z[i] = sqrt(nf_x[i]*nf_x[i] + nf_y[i]*nf_y[i]);
		        }

		        dA[i] = dr/2.;

      			// Calc the center of the element
      			ibm->cent_x[i]= (x_bp[n1e]+x_bp[n2e]+x_bp[n3e])/3.;
		        ibm->cent_y[i]= (y_bp[n1e]+y_bp[n2e]+y_bp[n3e])/3.;
		        ibm->cent_z[i]= (z_bp[n1e]+z_bp[n2e]+z_bp[n3e])/3.;
		}


    		ibm->nf_x = nf_x; ibm->nf_y = nf_y;  ibm->nf_z = nf_z;

    		ibm->dA = dA;
    		ibm->nt_x = nt_x; ibm->nt_y = nt_y;  ibm->nt_z = nt_z;
	        ibm->ns_x = ns_x; ibm->ns_y = ns_y;  ibm->ns_z = ns_z;

	        MPI_Bcast(ibm->nv1, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nv2, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nv3, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);

    	        MPI_Bcast(ibm->nf_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nf_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nf_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

    		MPI_Bcast(ibm->dA, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

    		MPI_Bcast(ibm->nt_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->nt_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->nt_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

    		MPI_Bcast(ibm->ns_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->ns_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->ns_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

    		MPI_Bcast(ibm->cent_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->cent_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->cent_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
    		int ti=0;
    		  	}
  	else if (rank) {
    		MPI_Bcast(&(n_v), 1, MPI_INT, 0, PETSC_COMM_WORLD);
	        ibm->n_v = n_v;

    	       	MPI_Bcast(&(ibm->n_elmt), 1, MPI_INT, 0, PETSC_COMM_WORLD);
	       	n_elmt = ibm->n_elmt;

	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp));

		x_bp = ibm->x_bp;	// added by seokkoo 03.04.2009
		y_bp = ibm->y_bp;
		z_bp = ibm->z_bp;

	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_i));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_i));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_i));

	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_o));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_o));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_o));

	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp0));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp0));
	        PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp0));


	        PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->u));
	        PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->uold));
	        PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->urm1));

		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->tmprt));


		PetscMalloc(n_v*sizeof(int), &ibm->node_boundary); // xiaolei add
      		PetscMalloc(n_elmt*sizeof(int), &ibm->elmt_boundary); // xiaolei add

     	 	PetscMalloc(3*n_elmt*sizeof(PetscReal), &(ibm->error_force)); // xiaolei add


		// fsi_fem
      		PetscMalloc(n_v*sizeof(Cmpnts), &(ibm->u_pre));

      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->x_bp_pre));
      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->y_bp_pre));
      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->z_bp_pre));


		// Error Vector
		PetscMalloc(9*n_v*sizeof(PetscReal), &(ibm->error));

	      	PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_x);
      		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_y);
      		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_z);

      		// xiaolei add
      		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_x_pre);
      		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_y_pre);
      		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->Traction_z_pre);

      		// xiaolei added
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Tx);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Ty);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumP_Tz);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->AbsP_tau);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Tx);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Ty);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->SumM_Tz);
		PetscMalloc(n_elmt*sizeof(PetscReal), &ibm->AbsM_tau);

	        for (i=0; i<n_v; i++) {

			ibm->tmprt[i] = 0.;

	  		ibm->u[i].x = 0.;
  		        ibm->u[i].y = 0.;
	    	        ibm->u[i].z = 0.;

		        ibm->uold[i].x = 0.;
		        ibm->uold[i].y = 0.;
	   	        ibm->uold[i].z = 0.;

	  	        ibm->urm1[i].x = 0.;
		        ibm->urm1[i].y = 0.;
	 	        ibm->urm1[i].z = 0.;
	        }

	        MPI_Bcast(ibm->x_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->y_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->z_bp_i, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->x_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->y_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->z_bp0, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->x_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->y_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->z_bp, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);

 	        MPI_Bcast(ibm->x_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	       	MPI_Bcast(ibm->y_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);
	       	MPI_Bcast(ibm->z_bp_o, n_v, MPIU_REAL, 0, PETSC_COMM_WORLD);



	       	PetscMalloc(n_elmt*sizeof(int), &nv1);
	       	PetscMalloc(n_elmt*sizeof(int), &nv2);
	       	PetscMalloc(n_elmt*sizeof(int), &nv3);

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nf_x);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nf_y);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nf_z);

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &dA);

   	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nt_x);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nt_y);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &nt_z);

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &ns_x);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &ns_y);
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &ns_z);

	       	ibm->nv1 = nv1; ibm->nv2 = nv2; ibm->nv3 = nv3;
	       	ibm->nf_x = nf_x; ibm->nf_y = nf_y; ibm->nf_z = nf_z;

	       	ibm->dA = dA;
	       	ibm->nt_x = nt_x; ibm->nt_y = nt_y;  ibm->nt_z = nt_z;
	       	ibm->ns_x = ns_x; ibm->ns_y = ns_y;  ibm->ns_z = ns_z;

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_x));
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_y));
	       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->cent_z));

      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_x));
      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_y));
      		PetscMalloc(n_v*sizeof(PetscReal), &(ibm->disp_z));
      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->disp_theta));

	       	// rotor model
	       	//


        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Urelmag_mean));

        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Fr_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ft_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Fa_mean));


        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->AOAAOA_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->FFt_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->FFa_mean));



        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ur_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ut_mean));
        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ua_mean));

        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->AOA_mean));



        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Urelmag));

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_x));
      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_y));
      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->F_lagr_z));

	       	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_x));
      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_y));
      		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dF_lagr_z));

	        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_x));
	        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_y));
	        PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->U_lagr_z));


        	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Ftmprt_lagr));
		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->Tmprt_lagr));




	        PetscMalloc(n_elmt*sizeof(int), &(ibm->i_min));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->j_min));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->k_min));

	        PetscMalloc(n_elmt*sizeof(int), &(ibm->i_max));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->j_max));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->k_max));

		PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhx));
          	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhy));
              	PetscMalloc(n_elmt*sizeof(PetscReal), &(ibm->dhz));

        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->angle_attack));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->angle_twist));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->chord_blade));

                PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CD));
                PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CL));

	        PetscMalloc(n_elmt*sizeof(int), &(ibm->iIP_min));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->jIP_min));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->kIP_min));

	        PetscMalloc(n_elmt*sizeof(int), &(ibm->iIP_max));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->jIP_max));
	        PetscMalloc(n_elmt*sizeof(int), &(ibm->kIP_max));

	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->random_color));

	       	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_x));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_y));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->centIP_z));

	       	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_x));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_y));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->U_IPlagr_z));


        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Nut_lagr));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->dh_IP));

	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->P_IPlagr));
	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_x));
	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_y));
	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->Shear_lagr_z));

	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_x));
	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_y));
	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->ShearDesired_lagr_z));

        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_x));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_y));
        	PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->UU_lagr_z));

	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Urel));
	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Uinduced));
	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->circulation));
	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->liftdirection));
	        PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->rotationdirection));

	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Urel_mean));
	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->Uinduced_mean));
	      	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->circulation_mean));
	       	PetscMalloc(ibm->n_elmt*sizeof(Cmpnts), &(ibm->liftdirection_mean));


	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->circulationSpecified));




	        PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->frictionfactor));

		PetscMalloc(n_elmt*sizeof(int), &(ibm->color));

		PetscMalloc(n_elmt*sizeof(int), &(ibm->s2l));

		PetscMalloc(ibm->n_elmt*sizeof(PetscReal), &(ibm->CTdisk));

	    	MPI_Bcast(ibm->nv1, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);
    		MPI_Bcast(ibm->nv2, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nv3, n_elmt, MPI_INT, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->nf_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nf_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nf_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->dA, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->nt_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nt_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->nt_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->ns_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->ns_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->ns_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);

	        MPI_Bcast(ibm->cent_x, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->cent_y, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	        MPI_Bcast(ibm->cent_z, n_elmt, MPIU_REAL, 0, PETSC_COMM_WORLD);
	}
	PetscPrintf(PETSC_COMM_WORLD, "Read %s !\n", fname);

    	for (i=0; i<ibm->n_elmt; i++) {
		ibm->Urelmag_mean[i]=0.0;
		ibm->Fr_mean[i]=0.0;
		ibm->Ft_mean[i]=0.0;
		ibm->Fa_mean[i]=0.0;
		ibm->Ur_mean[i]=0.0;
		ibm->Ut_mean[i]=0.0;
		ibm->Ua_mean[i]=0.0;
		ibm->AOA_mean[i]=0.0;
		ibm->color[i]=0;

		ibm->FFa_mean[i]=0.0;
		ibm->FFt_mean[i]=0.0;
		ibm->AOAAOA_mean[i]=0.0;

		ibm->Urel_mean[i].x=0.0;
		ibm->Urel_mean[i].y=0.0;
		ibm->Urel_mean[i].z=0.0;


		ibm->Uinduced_mean[i].x=0.0;
		ibm->Uinduced_mean[i].y=0.0;
		ibm->Uinduced_mean[i].z=0.0;

		ibm->circulation_mean[i].x=0.0;
		ibm->circulation_mean[i].y=0.0;
		ibm->circulation_mean[i].z=0.0;


		ibm->liftdirection_mean[i].x=0.0;
		ibm->liftdirection_mean[i].y=0.0;
		ibm->liftdirection_mean[i].z=0.0;


		ibm->disp_theta[i]=0;

	}

    	for (i=0; i<ibm->n_v; i++) {
		ibm->disp_x[i]=0;
		ibm->disp_y[i]=0;
		ibm->disp_z[i]=0;
	}


    	for (i=0; i<ibm->n_elmt; i++)
	{
		ibm->Traction_x[i] = 0;
		ibm->Traction_y[i] = 0;
		ibm->Traction_z[i] = 0;

		ibm->Traction_x_pre[i] = 0;
		ibm->Traction_y_pre[i] = 0;
		ibm->Traction_z_pre[i] = 0;
	}

	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: Surface_read_ucd  \n");

  	return(0);


}
//add over

PetscErrorCode velocity_initialize(IBMNodes *ibm){
    int rank,i;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
    for (i=0; i<ibm->n_elmt; i++) {
        ibm->U_lagr_x[i] = 0.001;
        ibm->U_lagr_y[i] = 0.001;
        ibm->U_lagr_z[i] = 0.5;
    }
    PetscPrintf(PETSC_COMM_WORLD, "------ velocity_initialize");
}


PetscErrorCode Force_initialize(IBMNodes *ibm){
    int rank,i;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
    for (i=0; i<ibm->n_elmt; i++) {
        ibm->F_lagr_x[i] = 2.;
        ibm->F_lagr_y[i] = 2.;
        ibm->F_lagr_z[i] = 10.;
    }
    PetscPrintf(PETSC_COMM_WORLD, "------ Force_initialize");
}


//wzw add : map the velocity from Eulerian grid to WriteData
PetscErrorCode Set_Vertices(PreciceInfo *precice,IBMNodes *ibm){
	int rank,i;
  	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
	PetscInt dim = precice->dim;
	PetscInt n_elmt = ibm->n_elmt;

  	precice->numberOfVertices = n_elmt;//the number of face centers
		
	PetscMalloc(precice->numberOfVertices * dim * sizeof(PetscReal), &(precice->vertices));
	PetscMalloc(precice->numberOfVertices * dim * sizeof(PetscReal), &(precice->readData));
	PetscMalloc(precice->numberOfVertices * dim * sizeof(PetscReal), &(precice->writeData));
	PetscMalloc(precice->numberOfVertices * sizeof(PetscInt), &(precice->vertexIDs));

	for (i=0; i<precice->numberOfVertices; i++) {
			precice->vertices[i*dim] = ibm->cent_x[i];
			precice->vertices[i*dim+1] = ibm->cent_y[i];
			precice->vertices[i*dim+2] = ibm->cent_z[i];

			precice->readData[i*dim] = 0.0;
			precice->readData[i*dim+1] = 0.0;
			precice->readData[i*dim+2] = 0.0;

			precice->writeData[i*dim] = 0.0;
			precice->writeData[i*dim+1] = 0.0;
			precice->writeData[i*dim+2] = 0.0;
		}
	
	precicec_setMeshVertices(precice->meshID, precice->numberOfVertices, precice->vertices,
                               precice->vertexIDs);

	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: Set_Vertices to initialize vertices, readData and writeData in precice \n");

}


//wzw add : map the velocity from Eulerian grid to WriteData
PetscErrorCode Map_WriteData(PreciceInfo *precice, IBMNodes *ibm) {
	int rank, i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
	PetscInt dim = precice->dim;

		for (i = 0; i < precice->numberOfVertices; i++) {
			precice->writeData[i * dim] = ibm->U_lagr_x[i];
			precice->writeData[i * dim + 1] = ibm->U_lagr_y[i];
			precice->writeData[i * dim + 2] = ibm->U_lagr_z[i];
		}
		
	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: Map_WriteData   \n");
}


//wzw add : map the Force from ReadData to ibm
PetscErrorCode Map_ReadData(PreciceInfo *precice, IBMNodes *ibm){
	int rank,i;
  	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
	PetscInt dim = precice->dim;

	for (i=0; i<precice->numberOfVertices; i++) {
			ibm->F_lagr_x[i] = precice->readData[i*dim];
			ibm->F_lagr_y[i] = precice->readData[i*dim+1];
			ibm->F_lagr_z[i] = precice->readData[i*dim+2];
	}

	PetscPrintf(PETSC_COMM_WORLD," ----------------Call: Map_ReadData   \n");
}


//wzw add for read surface data Preprocess
PetscErrorCode Interfaces_Preprocess(PetscReal NumberOfSurfaces,PetscReal cl,IBMNodes *VFS_OutPut, FSInfo *FSI_OutPut, double reflength_wt){  
	int my_rank,i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &my_rank);
	PetscInt ibi; // Add the missing variable declaration for ibi

    char path[256];
    sprintf(path, ".");

	if (!my_rank) {
		FILE *fd;
		char str[256];
		sprintf(str, "%s/Interfaces.inp", path);//store the interface center crood
		fd = fopen(str, "r");
		if(!fd) PetscPrintf(PETSC_COMM_WORLD, "Error: cannot open %s !\n", str),exit(0);
		for (ibi=0;ibi<NumberOfSurfaces;ibi++) {
			fscanf(fd, "%le %le %le", &(FSI_OutPut[ibi].x_c), &(FSI_OutPut[ibi].y_c), &(FSI_OutPut[ibi].z_c));
			FSI_OutPut[ibi].x_c=FSI_OutPut[ibi].x_c/cl;
			FSI_OutPut[ibi].y_c=FSI_OutPut[ibi].y_c/cl;
			FSI_OutPut[ibi].z_c=FSI_OutPut[ibi].z_c/cl;
			MPI_Bcast(&(FSI_OutPut[ibi].x_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
			MPI_Bcast(&(FSI_OutPut[ibi].y_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
			MPI_Bcast(&(FSI_OutPut[ibi].z_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
			PetscPrintf(PETSC_COMM_WORLD, "The normal center for %d th Surface %f %f %f \n", ibi, (FSI_OutPut[ibi].x_c), (FSI_OutPut[ibi].y_c), (FSI_OutPut[ibi].z_c));
		}
		for (ibi=0;ibi<NumberOfSurfaces;ibi++) {
			FSI_OutPut[ibi].x_c0=FSI_OutPut[ibi].x_c;
			FSI_OutPut[ibi].y_c0=FSI_OutPut[ibi].y_c;
			FSI_OutPut[ibi].z_c0=FSI_OutPut[ibi].z_c;
		}
		fclose(fd);
	}
	else {
		for (ibi=0;ibi<NumberOfSurfaces;ibi++) {
			MPI_Bcast(&(FSI_OutPut[ibi].x_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
			MPI_Bcast(&(FSI_OutPut[ibi].y_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
			MPI_Bcast(&(FSI_OutPut[ibi].z_c), 1, MPIU_REAL, 0, PETSC_COMM_WORLD);
		}
		for (ibi=0;ibi<NumberOfSurfaces;ibi++) {
			FSI_OutPut[ibi].x_c0=FSI_OutPut[ibi].x_c;
			FSI_OutPut[ibi].y_c0=FSI_OutPut[ibi].y_c;
			FSI_OutPut[ibi].z_c0=FSI_OutPut[ibi].z_c;
		}
	}

	double reflength = reflength_wt;
	for (i=0;i<NumberOfSurfaces;i++){
		char fname[80];
		sprintf(fname,"InterfacesUCD-%d",i+1);
		Surface_read_ucd(&(VFS_OutPut[i]), i, &(FSI_OutPut[i]),fname, reflength);

	}
	PetscPrintf(PETSC_COMM_WORLD, "Interfaces pre process over!\n");

}
//add over

//wzw add ： Initialize the data and file for precice coupling
PetscErrorCode Initialize_PRECICE(PreciceInfo *precice_list,PetscReal NumberOfSurfaces,PetscReal cl,IBMNodes *VFS_OutPut, FSInfo *FSI_OutPut, double reflength_wt,PetscInt ti){
	
    PetscBool flag;

	int Process_rank,Process_size,i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &Process_rank);
	MPI_Comm_size(PETSC_COMM_WORLD, &Process_size);
	PetscPrintf(PETSC_COMM_WORLD, "-----Initialize precice is coupling\n");

	PetscInt interfaces;//the number of interfaces in coupling,in my case is 1
	interfaces = NumberOfSurfaces;


	PetscInt dimensions;//the dimension of the case, in my case is 3
	dimensions = precicec_getDimensions();

	PetscPrintf(PETSC_COMM_WORLD, "-----Coupling dimensions is equal to :%d\n", dimensions);
	
	Interfaces_Preprocess(NumberOfSurfaces,cl,VFS_OutPut,FSI_OutPut,reflength_wt);

	PetscPrintf(PETSC_COMM_WORLD, "-----Reading Interface files is over!!!\n");


	for(i = 1;i<interfaces+1;i++){
		int index = i-1;
		precice_list[index].dim = dimensions;

		PetscPrintf(PETSC_COMM_WORLD, "-----show input info :  Index = index-%d----------\n", index);

		char filen[256];
		sprintf(filen, "Interface%d.dat", i);
		PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, filen, PETSC_TRUE);

		char meshName[256]; // Initialize the meshName variable
		PetscOptionsGetString(NULL,PETSC_NULL, "-meshName", meshName, 256, PETSC_NULL); // Assign the value to meshName
		precice_list[index].meshName = strdup(meshName); // Assign the value of meshName to precice_list[index].meshName
		PetscPrintf(PETSC_COMM_WORLD, "------meshName =  %s !!!--------\n", meshName);


		char dataName[256]; // Initialize the dataName variable
		PetscOptionsGetString(NULL,PETSC_NULL, "-writeData", dataName, 256, PETSC_NULL); // Assign the value to dataName
		precice_list[index].writeDataName = strdup(dataName); // Assign the value of dataName to precice_list[index].dataName
		PetscPrintf(PETSC_COMM_WORLD, "------writeData =  %s !!!--------\n", dataName);

		char readDataName[256]; // Initialize the readDataName variable
		PetscOptionsGetString(NULL,PETSC_NULL, "-readData", readDataName, 256, PETSC_NULL); // Assign the value to readDataName
		precice_list[index].readDataName = strdup(readDataName); // Assign the value of readDataName to precice_list[index].readDataName
		PetscPrintf(PETSC_COMM_WORLD, "------readData =  %s !!!--------\n", readDataName);

		char loctions[256]; // Initialize the readDataName variable
		PetscOptionsGetString(NULL,PETSC_NULL, "-locations", loctions, 256, PETSC_NULL); // Assign the value to readDataName
		PetscPrintf(PETSC_COMM_WORLD, "------locations =  %s !!!--------\n", loctions);

		precice_list[index].location = strdup(loctions); // Assign the value of readDataName to precice_list[index].readDataName
		PetscPrintf(PETSC_COMM_WORLD, "--------Read Interface%d.dat is over !!!!!\n", i);

		index++;


				/*
		locations = "faceCenters" or "faceNodes"
		DataName = "Displacement" or "Velocity" or "Pressure" or "Stress" or "Force" or "None", None means no data
		*/
	}


	for(i = 1;i<interfaces+1;i++){
		int index = i-1;
		
		PetscPrintf(PETSC_COMM_WORLD, "----- ======= ======= ======= ====== ======= ======= --------\n");
		PetscPrintf(PETSC_COMM_WORLD, "-----show precice meshiD info : Index = index-%d----------\n", index);


		precice_list[index].meshID = precicec_getMeshID(precice_list[index].meshName);

		PetscPrintf(PETSC_COMM_WORLD, "-----Coupling Interface%d----------\n", i);
		PetscPrintf(PETSC_COMM_WORLD, "-----Coupling meshName is equal to : %s !!!--------\n", precice_list[index].meshName);
		PetscPrintf(PETSC_COMM_WORLD, "-----Coupling readDataName is equal to : %s !!!--------\n", precice_list[index].readDataName);
		PetscPrintf(PETSC_COMM_WORLD, "-----Coupling writeDataName is equal to : %s !!!--------\n", precice_list[index].writeDataName);
		PetscPrintf(PETSC_COMM_WORLD, "-----Coupling meshID is equal to :%d !!!--------\n", precice_list[index].meshID);

		PetscPrintf(PETSC_COMM_WORLD, "----- ======= ======= ======= ====== ======= ======= --------\n");

		Set_Vertices(&(precice_list[index]),&(VFS_OutPut[index]));
	}

	for(i = 1;i<interfaces+1;i++){
		int index = i-1;
		PetscPrintf(PETSC_COMM_WORLD, "----- ======= ======= ======= ====== ======= ======= --------\n");
		PetscPrintf(PETSC_COMM_WORLD, "-----show precice dataID info : = index-%d----------\n", index);

		if(strcmp(precice_list[index].readDataName, "None") ==0){
			PetscPrintf(PETSC_COMM_WORLD, "-----There is no readDataName  !!!--------\n");
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling readDataName is equal to : %s !!!--------\n", precice_list[index].readDataName);
		}
		else{
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling readDataName is equal to : %s !!!--------\n", precice_list[index].readDataName);
			precice_list[index].readDataID = precicec_getDataID(precice_list[index].readDataName, precice_list[index].meshID);
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling readDataID is equal to : %d !!!--------\n", precice_list[index].readDataID);
		} 


		if(strcmp(precice_list[index].writeDataName, "None")==0){
			PetscPrintf(PETSC_COMM_WORLD, "-----There is no writeDataName  !!!--------\n");
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling writeDataName is equal to : %s !!!--------\n", precice_list[index].writeDataName);

		}
		else{
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling writeDataName is equal to : %s !!!--------\n", precice_list[index].writeDataName);
			precice_list[index].writeDataID = precicec_getDataID(precice_list[index].writeDataName, precice_list[index].meshID);
			PetscPrintf(PETSC_COMM_WORLD, "-----Coupling writeDataID is equal to : %d !!!--------\n", precice_list[index].writeDataID);

		}
		PetscPrintf(PETSC_COMM_WORLD, "-----Interface%d\n meshName:%s\n mesh Initialized!!!\n", i, precice_list[index].meshName);
		
	}
	
	PetscPrintf(PETSC_COMM_WORLD, "-----Call Initialize_PRECICE() !!! ----------------\n");
}

//wzw add : cal the data and output file for precice coupling
PetscErrorCode Read_Process_PRECICE(PreciceInfo *precice_list,PetscReal NumberOfSurfaces,IBMNodes *VFS_OutPut ,PetscInt ti,PetscReal dt){
	
	int Process_rank,Process_size,i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &Process_rank);
	MPI_Comm_size(PETSC_COMM_WORLD, &Process_size);
	PetscPrintf(PETSC_COMM_WORLD, "-----Read Process precice is coupling\n");

	PetscInt interfaces;//the number of interfaces in coupling,in my case is 1
	interfaces = NumberOfSurfaces;

	if(precicec_isReadDataAvailable()){

		PetscPrintf(PETSC_COMM_WORLD, "-----Read Data is available!!!\n");

		for(i = 0;i<interfaces;i++){
			
			if(strcmp(precice_list[i].readDataName, "None") != 0){
				PetscPrintf(PETSC_COMM_WORLD, "-----Interface-%d\n is reading data\n", i+1);
				precicec_readBlockVectorData(precice_list[i].readDataID, precice_list[i].numberOfVertices, precice_list[i].vertexIDs, precice_list[i].readData);
				Map_ReadData(&(precice_list[i]),&(VFS_OutPut[i]));
				//Force_surface_out(&(VFS_OutPut[i]), ti, i);
			}
			else{
				PetscPrintf(PETSC_COMM_WORLD, "-----Interface-%d\n not read data\n", i+1);
			}

		}

	}
	else{
		PetscPrintf(PETSC_COMM_WORLD, "-----Read Data is not available!!!\n");
	}

	for(i = 0;i<interfaces;i++){
		Force_surface_out(&(VFS_OutPut[i]), ti, i);
	} 

}


//wzw add : cal the data and output file for precice coupling
PetscErrorCode Write_Process_PRECICE(PreciceInfo *precice_list,PetscReal NumberOfSurfaces,IBMNodes *VFS_OutPut ,PetscInt ti,PetscReal dt){
	
	int Process_rank,Process_size,i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &Process_rank);
	MPI_Comm_size(PETSC_COMM_WORLD, &Process_size);
	PetscPrintf(PETSC_COMM_WORLD, "-----Write Process precice is coupling\n");

	PetscInt interfaces;//the number of interfaces in coupling,in my case is 1
	interfaces = NumberOfSurfaces;

	if(precicec_isWriteDataRequired(dt)){

		PetscPrintf(PETSC_COMM_WORLD, "-----Write Data is available!!!\n");

		for(i = 0;i<interfaces;i++){
			
			if(strcmp(precice_list[i].writeDataName, "None") != 0){
				PetscPrintf(PETSC_COMM_WORLD, "-----Interface-%d\n is Writing data\n", i+1);
                velocity_initialize(&(VFS_OutPut[i]));
				Map_WriteData(&(precice_list[i]),&(VFS_OutPut[i]));
				precicec_writeBlockVectorData(precice_list[i].writeDataID, precice_list[i].numberOfVertices, precice_list[i].vertexIDs, precice_list[i].writeData);
			}
			else{
				PetscPrintf(PETSC_COMM_WORLD, "-----Interface-%d\n not Write data\n", i+1);
			}
		}
	}
	else{
		PetscPrintf(PETSC_COMM_WORLD, "-----write Data is not available!!!\n");
	}

	for(i = 0;i<interfaces;i++){
		velocity_surface_out(&(VFS_OutPut[i]), ti, i);
	} 
	//end of write data!!!

}

PetscErrorCode Write_DataFile(PreciceInfo *precice_list,PetscReal NumberOfSurfaces,IBMNodes *VFS_OutPut ,PetscInt ti){
	int Process_rank,Process_size,i;
	MPI_Comm_rank(PETSC_COMM_WORLD, &Process_rank);
	MPI_Comm_size(PETSC_COMM_WORLD, &Process_size);
	PetscPrintf(PETSC_COMM_WORLD, "-----Write Process precice is coupling\n");

	PetscInt interfaces;//the number of interfaces in coupling,in my case is 1
	interfaces = NumberOfSurfaces;

	for(i = 0;i<interfaces;i++){
		PetscPrintf(PETSC_COMM_WORLD, "-----Interface-%d\n is Writing data\n", i+1);
        velocity_initialize(&(VFS_OutPut[i]));
		Force_initialize(&(VFS_OutPut[i]));
		velocity_surface_out(&(VFS_OutPut[i]), ti, i);
		Force_surface_out(&(VFS_OutPut[i]), ti, i);
	} 
	//end of write data!!!
}


int main(int argc, char **argv)
{
    PetscErrorCode ierr;

    ierr = PetscInitialize(&argc, &argv, (char *)0, NULL);
    if (ierr != 0) {
        PetscPrintf(PETSC_COMM_WORLD, "PetscInitialize failed!\n");
        return ierr;
    }

    PetscReal dt = 0.0;
    PetscInt ti,tistart=0;
    PetscInt dimensions = -1;

    /*WZW: Coupling1.0 2024/1/19*/
    char ConfigFileName[256];
    char ParticipantName[256];
    PetscInt CoupleOtherSolver=0;
    PetscInt  NumInterFaces = 0;//VFS-WIND to other solver(output Velocity and Pressure), or receive Forces from other solver
    PetscReal precice_dt = 1000.0;//dt for precice
    PetscReal solver_dt = 10000.0;
    PetscReal cl_center = 1.;
	PetscReal cl_nodes = 1.;

    //WZW: add for coupling
	IBMNodes 	  *VFS_OutPut;//interfaces in VFS
	FSInfo		  *FSI_OutPut;//interfaces in VFS
	PreciceInfo *precice_list;//precice in VFS



    PetscOptionsInsertFile(PETSC_COMM_WORLD, NULL, "control.dat", PETSC_TRUE);
    //WZW: add---start
	PetscOptionsGetInt(NULL,PETSC_NULL, "-isCoupling", &CoupleOtherSolver, PETSC_NULL);
	PetscOptionsGetInt(NULL,PETSC_NULL, "-num_interfaces", &NumInterFaces, PETSC_NULL);
	PetscOptionsGetReal(NULL, PETSC_NULL, "-dt", &solver_dt, PETSC_NULL);
	PetscOptionsGetString(NULL, PETSC_NULL, "-preciceConfig", ConfigFileName, 256, PETSC_NULL);
    PetscOptionsGetString(NULL, PETSC_NULL, "-participant", ParticipantName, 256, PETSC_NULL);
	//add---end

    // WZW: add for couling output
	if(CoupleOtherSolver>0){
		if (NumInterFaces>0){
			PetscMalloc(NumInterFaces*sizeof(IBMNodes), &VFS_OutPut);
			PetscMalloc(NumInterFaces*sizeof(FSInfo), &FSI_OutPut);
			PetscMalloc(NumInterFaces*sizeof(PreciceInfo), &precice_list);
			PetscPrintf(PETSC_COMM_WORLD, "---- Coupling is open, There are %d interfaces for VFS-Wind !\n", NumInterFaces);
		}
		else{
			PetscPrintf(PETSC_COMM_WORLD, "---- Coupling is open, But there is no interface for VFS-Wind !\n");
		}
	}
	else{
		PetscPrintf(PETSC_COMM_WORLD, "---- Coupling is closed !\n");
	}
	//end

    PetscMPIInt commRank = -1, commSize = -1;
    MPI_Comm_rank(PETSC_COMM_WORLD, &commRank);
    MPI_Comm_size(PETSC_COMM_WORLD, &commSize);

    if (CoupleOtherSolver>0){
		ti = 0;

		PetscOptionsGetReal(NULL,PETSC_NULL, "-cl_center_Precice", &cl_center, PETSC_NULL);
		PetscOptionsGetReal(NULL,PETSC_NULL, "-cl_nodes_Precice", &cl_nodes, PETSC_NULL);
		PetscPrintf(PETSC_COMM_WORLD, "-----DEBUG 1 is ok!!! ------\n");

		precicec_createSolverInterface(ParticipantName, ConfigFileName, commRank, commSize);
		PetscPrintf(PETSC_COMM_WORLD, "-----DEBUG 2 is ok!!! ------\n");

		Initialize_PRECICE(precice_list,NumInterFaces,cl_center,VFS_OutPut, FSI_OutPut,cl_nodes,ti);
		PetscPrintf(PETSC_COMM_WORLD, "-----Initialize precice coupling is ok!!! ------\n");

		Write_DataFile(precice_list,NumInterFaces,VFS_OutPut ,ti);
		PetscPrintf(PETSC_COMM_WORLD, "-----DATA file write over!!!! \n");

		precice_dt = precicec_initialize();
		PetscPrintf(PETSC_COMM_WORLD, "Precice precice_dt is equl to %le\n", precice_dt);
		PetscBarrier(PETSC_NULL);
	}


    while (precicec_isCouplingOngoing()) {

        if(CoupleOtherSolver>0){
				PetscReal min_dt;
				if(precice_dt>solver_dt){
					min_dt = solver_dt;
				}
				else{
					min_dt = precice_dt;
				}
                				dt = min_dt;
				PetscPrintf(PETSC_COMM_WORLD, "Precice min dt is equl to %le\n", dt);
				Read_Process_PRECICE(precice_list,NumInterFaces,VFS_OutPut ,ti,dt);
			}
			//end
			/*wzw add for precice write data*/
		if(CoupleOtherSolver>0){
				Write_Process_PRECICE(precice_list,NumInterFaces,VFS_OutPut ,ti,dt);
				precice_dt = precicec_advance(dt);
			}
			//end

        ti++;
    }

    precicec_finalize();
    PetscFree(VFS_OutPut);
    PetscFree(FSI_OutPut);
    PetscFree(precice_list);

    ierr = PetscFinalize();
    if (ierr != 0) {
        PetscPrintf(PETSC_COMM_WORLD, "PetscFinalize failed!\n");
        return ierr;
    }

    return 0;
}