
// set_parity_conditions_from_symbolic_dot_products():
// Evaluate dot products needed for setting parity
//     conditions at a given point (xx0,xx1,xx2),
//     using C code generated by NRPy+ function
//     parity_conditions_symbolic_dot_products().
void eval_symbolic_dot_products_to_set_parity_conditions(const paramstruct *restrict params,  REAL parity[10],
                                    const REAL xx0,const REAL xx1,const REAL xx2,
                                    const REAL xx0_inbounds,const REAL xx1_inbounds,const REAL xx2_inbounds) {
#include "RELATIVE_PATH__set_Cparameters.h" /* Header file containing correct #include for set_Cparameters.h;
                                             * accounting for the relative path */
#include "parity_conditions_symbolic_dot_products.h"
}

void set_up__bc_gz_map_and_parity_condns(const paramstruct *restrict params,
                                         REAL *xx[3], gz_map *bc_gz_map,parity_condition *bc_parity_conditions) {
#include "RELATIVE_PATH__set_Cparameters.h" /* Header file containing correct #include for set_Cparameters.h;
                                             * accounting for the relative path */
  // xx[0][j] = xxmin[0] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*dxx0;
  // -> xxmin[0] = xx[0][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx0
  const REAL xxmin[3] = { xx[0][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx0,
                          xx[1][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx1,
                          xx[2][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx2 };
  //fprintf(stderr,"hey inside setbc: %e %e %e | %e %e\n",xxmin[0],xxmin[1],xxmin[2],xx[0][0],dxx0);
  LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,Nxx_plus_2NGHOSTS1,0,Nxx_plus_2NGHOSTS2) {
    // Step 1: Convert the (curvilinear) coordinate (x0,x1,x2) to Cartesian coordinates
    REAL xCart[3];
    EigenCoord_xxCart(params, xx, i0,i1,i2, xCart);
    REAL Cartx = xCart[0];
    REAL Carty = xCart[1];
    REAL Cartz = xCart[2];

    // Step 2: Find the (i0_inbounds,i1_inbounds,i2_inbounds) corresponding to the above Cartesian coordinate.
    //   If (i0_inbounds,i1_inbounds,i2_inbounds) is in a ghost zone, then it must equal (i0,i1,i2), and
    //      the point is an outer boundary point.
    //   Otherwise (i0_inbounds,i1_inbounds,i2_inbounds) is in the grid interior, and data at (i0,i1,i2)
    //      must be replaced with data at (i0_inbounds,i1_inbounds,i2_inbounds), but multiplied by the
    //      appropriate parity condition (+/- 1).
    REAL Cart_to_xx0_inbounds,Cart_to_xx1_inbounds,Cart_to_xx2_inbounds;
#include "EigenCoord_Cart_to_xx.h"
    int i0_inbounds = (int)( (Cart_to_xx0_inbounds - xxmin[0] - (1.0/2.0)*dxx0 + ((REAL)NGHOSTS)*dxx0)/dxx0 + 0.5 );
    int i1_inbounds = (int)( (Cart_to_xx1_inbounds - xxmin[1] - (1.0/2.0)*dxx1 + ((REAL)NGHOSTS)*dxx1)/dxx1 + 0.5 );
    int i2_inbounds = (int)( (Cart_to_xx2_inbounds - xxmin[2] - (1.0/2.0)*dxx2 + ((REAL)NGHOSTS)*dxx2)/dxx2 + 0.5 );

    // Step 2.a: (Sanity/validation check) Convert the interior point
    //           x0(i0_inbounds),x1(i1_inbounds),x2(i2_inbounds) to Cartesian coordinates,
    //           make sure that the Cartesian coordinate matches the Cartesian coordinate of
    //           x0(i0),x1(i1),x2(i2). If not, error out!
    REAL xCart_orig[3]; for(int ii=0;ii<3;ii++) xCart_orig[ii] = xCart[ii];
    EigenCoord_xxCart(params, xx, i0_inbounds,i1_inbounds,i2_inbounds, xCart);

//fprintf(stderr,"Cartesian agreement: ( %.15e %.15e %.15e ) ?= ( %.15e %.15e %.15e )\n",
// (double)xCart_orig[0],(double)xCart_orig[1],(double)xCart_orig[2],
// (double)xCart[0],(double)xCart[1],(double)xCart[2]);

#define EPS_ABS 1e-8
    if(fabs( (double)(xCart_orig[0] - xCart[0]) ) > EPS_ABS ||
       fabs( (double)(xCart_orig[1] - xCart[1]) ) > EPS_ABS ||
       fabs( (double)(xCart_orig[2] - xCart[2]) ) > EPS_ABS) {
      fprintf(stderr,"Error. Cartesian disagreement: ( %.15e %.15e %.15e ) != ( %.15e %.15e %.15e )\n",
             (double)xCart_orig[0],(double)xCart_orig[1],(double)xCart_orig[2],
             (double)xCart[0],(double)xCart[1],(double)xCart[2]);
      exit(1);
    }

    // Step 3: Set bc_gz_map and bc_parity_conditions.
    if(i0_inbounds-i0 == 0 && i1_inbounds-i1 == 0 && i2_inbounds-i2 == 0) {
      // Step 3.a: Iff we are on an outer boundary point or in the grid
      //           interior, i0_inbounds==i0, i1_inbounds==i1, and
      //           i2_inbounds==i2, and inner boundary conditions do not
      //           apply: set bc_gz_map to -1, and parity=1.
      bc_gz_map[IDX3S(i0,i1,i2)].i0=-1;
      bc_gz_map[IDX3S(i0,i1,i2)].i1=-1;
      bc_gz_map[IDX3S(i0,i1,i2)].i2=-1;
      for(int which_parity=0; which_parity<10; which_parity++) {
        bc_parity_conditions[IDX3S(i0,i1,i2)].parity[which_parity] = 1;
      }
    } else {
      // Step 3.b: If we are on an *inner* boundary point:
      // 1. Set bc_gz_map at (i0,i1,i2) to the point
      //    in the interior to which this boundary
      //    point maps, and
      // 2. Perform the unit vector dot products
      //    necessary to set all 10 possible parity
      //    conditions, calling function
      //    set_parity_from_unit_vector_dot_product()
      bc_gz_map[IDX3S(i0,i1,i2)].i0=i0_inbounds;
      bc_gz_map[IDX3S(i0,i1,i2)].i1=i1_inbounds;
      bc_gz_map[IDX3S(i0,i1,i2)].i2=i2_inbounds;
      const REAL xx0 = xx[0][i0];
      const REAL xx1 = xx[1][i1];
      const REAL xx2 = xx[2][i2];
      const REAL xx0_inbounds = xx[0][i0_inbounds];
      const REAL xx1_inbounds = xx[1][i1_inbounds];
      const REAL xx2_inbounds = xx[2][i2_inbounds];
      REAL REAL_parity_array[10];
      eval_symbolic_dot_products_to_set_parity_conditions(params,  REAL_parity_array,  xx0,xx1,xx2,
                                                          xx0_inbounds,xx1_inbounds,xx2_inbounds);
      for(int whichparity=0;whichparity<10;whichparity++) {
          //printf("Good? Parity %d evaluated to %e\n",whichparity,(double)REAL_parity_array[whichparity]);
          // Perform sanity check on parity array output: should be +1 or -1 to within 8 significant digits:
          if( (REAL_parity_array[whichparity]  > 0 && fabs(REAL_parity_array[whichparity] - (+1)) > 1e-8) ||
              (REAL_parity_array[whichparity] <= 0 && fabs(REAL_parity_array[whichparity] - (-1)) > 1e-8) ) {
              fprintf(stderr,"Error at point (%d %d %d); (%e %e %e); maps to (%e %e %e).\n",
                      i0,i1,i2, xx0,xx1,xx2, xx0_inbounds,xx1_inbounds,xx2_inbounds);
              fprintf(stderr,"Parity evaluated to %e , which is not within 8 significant digits of +1 or -1.\n",
                      REAL_parity_array[whichparity]);
              exit(1);
          }
          if(REAL_parity_array[whichparity] < 0.0) bc_parity_conditions[IDX3S(i0,i1,i2)].parity[whichparity] = -1;
          if(REAL_parity_array[whichparity] > 0.0) bc_parity_conditions[IDX3S(i0,i1,i2)].parity[whichparity] = +1;
      }
    }
  }
}
