path: root/src/patch/patch.cc

// patch.cc -- describes a coordinate/grid patch
// $Id$

//
// patch::patch
// patch::~patch
// z_patch::z_patch
// x_patch::x_patch
// y_patch::y_patch
//
// patch::setup_mirror_sym_ghost_zone
// patch::setup_periodic_BC_ghost_zone
// patch::setup_interpatch_ghost_zone
// patch::setup_other_patch_frontier
// patch::setup_patch_frontier
// patch::border_on_edge
//
// patch::write_grid3d_gridfn_at_iwr
// patch::write_grid3d_gridfn__adat
// patch::radially_extrapolate_gridfn_vector
// *** template instantiations for  patch::radially_extrapolate_gridfn_vector
//

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

#include <vector.h>

#include <jt/stdc.h>
#include <jt/util.h>
#include <jt/util++.hh>
#include <jt/misc++.hh>

#include "policy.hh"
#include "fp.hh"
#include "misc.hh"
#include "coords.hh"
#include "gfa.hh"
#include "gridfn.hh"
#include "grid.hh"
#include "fd_grid.hh"
#include "patch.hh"
#include "patch_bndry.hh"

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This function constructs a  patch  object.
//
patch::patch(patch_system &my_patch_system_in, int patch_number_in,
	     const char *name_in, bool is_plus_in, char ctype_in,
	     local_coords::coords_set coords_set_rho_in,
				      coords_set_sigma_in,
				      coords_set_tau_in,
	     const grid_arrays::array_pars& grid_array_pars_in,
	     const grid::grid_pars& grid_pars_in,
	     int N_gridfns_in)

	: fd_grid(grid_array_pars_in, gride_pars_in, N_gridfns_in),

	  my_patch_system_(my_patch_system_in), patch_number_(patch_number_in),
	  name_(name_in), is_plus_(is_plus_in), ctype_(ctype_in),

	  coords_set_rho_  (coords_set_rho_in  ),
	  coords_set_sigma_(coords_set_sigma_in),
	  coords_set_tau_  (coords_set_tau_in  ),

	  min_rho_patch_edge_  (*new patch_edge(*this, side_min, side_rho  )),
	  max_rho_patch_edge_  (*new patch_edge(*this, side_max, side_rho  )),
	  min_sigma_patch_edge_(*new patch_edge(*this, side_min, side_sigma)),
	  max_sigma_patch_edge_(*new patch_edge(*this, side_max, side_sigma)),

	  min_rho_patch_frontier_(NULL),
	  max_rho_patch_frontier_(NULL),
	  min_sigma_patch_frontier_(NULL),
	  max_sigma_patch_frontier_(NULL),

	  min_rho_ghost_zone_(NULL),
	  max_rho_ghost_zone_(NULL),
	  min_sigma_ghost_zone_(NULL),
	  max_sigma_ghost_zone_(NULL) // no comma

{ }

//*****************************************************************************

//
// This function destroys a  patch  object.
//
patch::~patch()
{
// no need to check for null pointers, since  delete NULL  is a silent no-op

delete max_sigma_ghost_zone_;
delete min_sigma_ghost_zone_;
delete max_rho_ghost_zone_;
delete min_rho_ghost_zone_;

delete max_sigma_patch_frontier_;
delete min_sigma_patch_frontier_;
delete max_rho_patch_frontier_;
delete min_rho_patch_frontier_;

delete max_sigma_patch_edge_;
delete min_sigma_patch_edge_;
delete max_rho_patch_edge_;
delete min_rho_patch_edge_;
}

//*****************************************************************************

//
// This function constructs a  z_patch  object.
//
z_patch::z_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in,
		 int N_gridfns_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'z',
		local_coords::set_mu, local_coords::set_nu,
		local_coords::set_phi,
		grid_array_pars_in, pars_in, N_gridfns_in)
{ }

//*****************************************************************************

//
// This function constructs an  x_patch  object.
//
x_patch::x_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in,
		 int N_gridfns_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'x',
		local_coords::set_nu, local_coords::set_phi,
		local_coords::set_mu,
		grid_array_pars_in, pars_in, N_gridfns_in)

//*****************************************************************************

//
// This function constructs a  y_patch  object.
//
y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
		 const char *name_in, bool is_plus_in,
		 const grid_arrays::array_pars& grid_array_pars_in,
		 const grid::grid_pars& grid_pars_in,
		 int N_gridfns_in)
	: patch(my_patch_system_in, patch_number_in,
		name_in, is_plus_in, 'y',
		local_coords::set_mu, local_coords::set_phi,
		local_coords::set_nu,
		grid_array_pars_in, pars_in, N_gridfns_in)

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This function sets up a mirror-symmetry angular border of this patch.
//
void patch::setup_mirror_sym_ghost_zone(const patch_edge& edge_in)
{
// make sure we belong to the right patch
assert(& edge_in.this_patch() == this);

symmetry_ghost_zone *temp = new symmetry_ghost_zone(edge_in);
symmetry_ghost_zones_array_[N_symmetry_ghost_zones_++] = temp;

// assignment converts to ptr to base class  ghost_zone::
// FIXME: the decode here is ugly; an array would be cleaner...
if (edge_in.edge_is_min())
   then if (edge_in.edge_is_rho())
	   then min_rho_ghost_zone_   = temp;
	   else min_sigma_ghost_zone_ = temp;
   else if (edge_in.edge_is_rho())
	   then max_rho_ghost_zone_   = temp;
	   else max_sigma_ghost_zone_ = temp;
}

//*****************************************************************************

//
// This function sets up a periodic-BCs angular border of this patch.
//
void patch::setup_periodic_BC_ghost_zone
	(const patch_edge& this_edge_in, const patch_edge& symmetry_edge_in,
	 int this_edge_sample_ipar, int symmetry_edge_sample_ipar,
	 bool ipar_map_sign)		// one of cpm_map::{plus,minus}_map
{
// make sure we belong to the right patch
assert(& this_edge_in.this_patch() == this);

symmetry_ghost_zone *temp
	= new symmetry_ghost_zone
		      (this_edge_in, symmetry_edge_in,
		       this_edge_sample_ipar, symmetry_edge_sample_ipar,
		       ipar_map_sign);
symmetry_ghost_zones_array_[N_symmetry_ghost_zones_++] = temp;

// assignment converts to ptr to base class  ghost_zone::
// FIXME: the decode here is ugly; an array would be cleaner...
if (this_edge_in.edge_is_min())
   then if (this_edge_in.edge_is_rho())
	   then min_rho_ghost_zone_   = temp;
	   else min_sigma_ghost_zone_ = temp;
   else if (this_edge_in.edge_is_rho())
	   then max_rho_ghost_zone_   = temp;
	   else max_sigma_ghost_zone_ = temp;
}

//*****************************************************************************

//
// This function sets up an interpatch angular border of this patch.
//
void patch::setup_interpatch_ghost_zone(const patch_edge& this_edge_in,
				    const patch_edge& other_edge_in)
{
// make sure we belong to the right patch
assert(& this_edge_in.this_patch() == this);

interpatch_ghost_zone *temp
	= new interpatch_ghost_zone(this_edge_in, other_edge_in);
interpatch_ghost_zones_array_[N_interpatch_ghost_zones_++] = temp;

// assignment converts to ptr to base class  ghost_zone::
// FIXME: the decode here is ugly; an array would be cleaner...
if (this_edge_in.edge_is_min())
   then if (this_edge_in.edge_is_rho())
	   then min_rho_ghost_zone_ = temp;
	   else min_sigma_ghost_zone_ = temp;
   else if (this_edge_in.edge_is_rho())
	   then max_rho_ghost_zone_ = temp;
	   else max_sigma_ghost_zone_ = temp;
}

//*****************************************************************************

//
// This function sets an angular frontier of *another* patch from
// which one of this patch's interpatch borders needs to interpolate
// data.
//
void patch::setup_other_patch_frontier
	(const patch_edge& this_edge_in,
	 const patch_edge& other_edge_in,
	 const interpolator_pars& interpolator_pars_in)
{
// make sure we belong to the right patch
assert(& this_edge_in.this_patch() == this);

const patch_frontier& temp
	= other_edge_in
	  .nonconst_this_patch()
	  .setup_patch_frontier(other_edge_in, interpolator_pars_in);

border_on_edge(this_edge_in)
.treat_me_as_nonconst_interpatch_ghost_zone()
.set_other_patch_frontier(temp);
}

//*****************************************************************************

//
// This function sets up an angular frontier of this patch.
//
const patch_frontier& patch::setup_patch_frontier
	(const patch_edge& edge_in,
	 const interpolator_pars& interpolator_pars_in)
{
// make sure we belong to the right patch
assert(& edge_in.this_patch() == this);

const patch_frontier *temp
	= new patch_frontier(edge_in, interpolator_pars_in);

// FIXME: the decode here is ugly; an array would be cleaner...
if (edge_in.edge_is_min())
   then if (edge_in.edge_is_rho())
	   then min_rho_patch_frontier_ = temp;
	   else min_sigma_patch_frontier_ = temp;
   else if (edge_in.edge_is_rho())
	   then max_rho_patch_frontier_ = temp;
	   else max_sigma_patch_frontier_ = temp;

return *temp;
}

//*****************************************************************************

//
// This function returns the ghost_zone object on a specified edge
// of this patch.
//
const ghost_zone& patch::border_on_edge(const patch_edge& edge_in)
	const
{
return edge_in.edge_is_min()
       ? ( edge_in.edge_is_rho() ? min_rho_ghost_zone()
				 : min_sigma_ghost_zone() )
       : ( edge_in.edge_is_rho() ? max_rho_ghost_zone()
				 : max_sigma_ghost_zone() );
}

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This function writes the values of a grid3d gridfn at a specified iwr
// to the already-open stdio stream  output_fp  (which is *not* closed
// when we're through), using the format
//	# dpx	dpy	gridfn
// Each set of  dpx = const lines is followed by a blank line.
// Each patch is followed by two blank lines.
//
void patch::write_grid3d_gridfn__adat(FILE *output_fp,
				      int gfn,
				      int iwr, bool want_borders)
	const
{
if (output_fp == NULL)
   then error_exit(PANIC_EXIT,
"***** patch::write_grid3d_gridfn__adat:\n"
"        got NULL output stream for gfn=%d iwr=%d!\n"
,
		   gfn, iwr);					/*NOTREACHED*/

fprintf(output_fp,
	"# %s patch: (dpx,dpy) = (%s,%s)\n",
	name(), defn_of_dpx(), defn_of_dpy());
fprintf(output_fp,
	"# dpx\tdpy\tgridfn\n");

	for (int irho = effective_min_irho(want_borders) ;
	     irho <= effective_max_irho(want_borders) ;
	     ++irho)
	{
		for (int isigma = effective_min_isigma(want_borders) ;
		     isigma <= effective_max_isigma(want_borders) ;
		     ++isigma)
		{
		fprintf(output_fp,
			"%g\t%g\t%.20e\n",
			double(dpx_of_irho_isigma(irho,isigma)),
			double(dpy_of_irho_isigma(irho,isigma)),
			double(gridfn__grid3d(gfn, iwr, irho, isigma)));
		}
	fprintf(output_fp, "\n");
	}

fprintf(output_fp, "\n");
}

//*****************************************************************************

//
// This function does radial extrapolation at a specified radial border,
// for all the gridfns in a specified gridfn vector.
// The loop nesting is cache-optimized for interleaved gridfns.
//
template <int pm>
void patch::radially_extrapolate_gridfn_vector(radial_extrapolation_type ret,
					       const gridfn_vector& gv)
{
	for (int irho = bordered_min_irho() ;
	     irho <= bordered_max_irho() ;
	     ++irho)
	{
	for (int isigma = bordered_min_isigma() ;
	     isigma <= bordered_max_isigma() ;
	     ++isigma)
	{
		for (int ii = 0 ; ii < gv.N_gridfns() ; ++ii)
		{
		const int gfn = gv.ith_gfn(ii);

		radially_extrapolate<pm>(ret, gfn, irho,isigma);
		}
	}
	}
}

//*****************************************************************************

//
// *** template instantiations for  patch::radially_extrapolate_gridfn_range
//

template
  void patch::radially_extrapolate_gridfn_vector<inner_boundary>
						(radial_extrapolation_type ret,
						 const gridfn_vector& gv);
template
  void patch::radially_extrapolate_gridfn_vector<outer_boundary>
						(radial_extrapolation_type ret,
						 const gridfn_vector& gv);