finish ghost-zone and patch-interp Jacobian support -- yea!!!!

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@598 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 94 insertions, 33 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index ecb397e..8516241 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -12,7 +12,9 @@
 // interpatch_ghost_zone::[min,max]_ipar
 // interpatch_ghost_zone::finish_setup
 // interpatch_ghost_zone::assert_fully_setup
+// interpatch_ghost_zone::verify_caller_wants_all_of_ghost_zone
 // interpatch_ghost_zone::synchronize
+// interpatch_ghost_zone::compute_Jacobian
 //
 
 #include <stdio.h>
@@ -123,7 +125,6 @@ void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 				      bool want_min_par_corner,
 				      bool want_non_corner,
 				      bool want_max_par_corner)
-	const
 {
 	for (int gfn = ghosted_min_gfn ; gfn <= ghosted_max_gfn ; ++gfn)
 	{
@@ -173,7 +174,8 @@ interpatch_ghost_zone::interpatch_ghost_zone(const patch_edge& my_edge_in,
 	  other_iperp_(NULL),
           min_ipar_used_(NULL),	max_ipar_used_(NULL),
 	  other_par_(NULL),
-	  interp_result_buffer_(NULL)
+	  interp_result_buffer_(NULL),
+	  Jacobian_oipar_posn_(NULL), Jacobian_buffer_(NULL) // no comma
 {
 //
 // verify that we have the expected relationships between
@@ -327,6 +329,8 @@ other_iperp_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 //
 interpatch_ghost_zone::~interpatch_ghost_zone()
 {
+delete Jacobian_buffer_;
+delete Jacobian_oipar_posn_;
 delete interp_result_buffer_;
 delete other_par_;
 delete max_ipar_used_;
@@ -387,10 +391,8 @@ return max_par_adjacent_ghost_zone().is_symmetry()
 void interpatch_ghost_zone::finish_setup(int interp_handle,
 					 int interp_par_table_handle)
 {
-const int other_min_iperp = std::min(other_iperp(min_iperp()),
-				     other_iperp(max_iperp()));
-const int other_max_iperp = std::max(other_iperp(min_iperp()),
-				     other_iperp(max_iperp()));
+min_other_iperp_ = std::min(other_iperp(min_iperp()), other_iperp(max_iperp()));
+max_other_iperp_ = std::max(other_iperp(min_iperp()), other_iperp(max_iperp()));
 
 
 //
@@ -400,16 +402,16 @@ const int other_max_iperp = std::max(other_iperp(min_iperp()),
 // ... we will pass these arrays by reference to the other patch's
 //     patch_interp object, with ipar being parindex there
 //
-int extreme_min_ipar = +INT_MAX;
-int extreme_max_ipar = -INT_MAX;
-min_ipar_used_ = new jtutil::array1d<int>(other_min_iperp, other_max_iperp);
-max_ipar_used_ = new jtutil::array1d<int>(other_min_iperp, other_max_iperp);
+extreme_min_ipar_ = INT_MAX;
+extreme_max_ipar_ = INT_MIN;
+min_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
+max_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
 	const int min_ipar_here = min_ipar(iperp);
 	const int max_ipar_here = max_ipar(iperp);
-	extreme_min_ipar = std::min(extreme_min_ipar, min_ipar_here);
-	extreme_max_ipar = std::max(extreme_max_ipar, max_ipar_here);
+	extreme_min_ipar_ = std::min(extreme_min_ipar_, min_ipar_here);
+	extreme_max_ipar_ = std::max(extreme_max_ipar_, max_ipar_here);
 	(*min_ipar_used_)(other_iperp(iperp)) = min_ipar_here;
 	(*max_ipar_used_)(other_iperp(iperp)) = max_ipar_here;
 	}
@@ -419,8 +421,8 @@ max_ipar_used_ = new jtutil::array1d<int>(other_min_iperp, other_max_iperp);
 // set up array giving other patch's par coordinate for interpolation
 //
 
-other_par_ = new jtutil::array2d<fp>(other_min_iperp, other_max_iperp,
-				     extreme_min_ipar, extreme_max_ipar);
+other_par_ = new jtutil::array2d<fp>(min_other_iperp_, max_other_iperp_,
+				     extreme_min_ipar_, extreme_max_ipar_);
 
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
@@ -449,16 +451,16 @@ other_par_ = new jtutil::array2d<fp>(other_min_iperp, other_max_iperp,
 // set up interpolation result buffer
 //
 interp_result_buffer_
-	= new jtutil::array3d<fp>
-		(my_patch().ghosted_min_gfn(), my_patch().ghosted_max_gfn(),
-		 other_min_iperp, other_max_iperp,
-		 extreme_min_ipar, extreme_max_ipar);
+	= new jtutil::array3d<fp>(my_patch().ghosted_min_gfn(),
+					  my_patch().ghosted_max_gfn(),
+				  min_other_iperp_, max_other_iperp_,
+				  extreme_min_ipar_, extreme_max_ipar_);
 
 //
 // construct the patch_interp object to interpolate from the *other* patch
 //
 other_patch_interp_ = new patch_interp(other_edge(),
-				       other_min_iperp, other_max_iperp,
+				       min_other_iperp_, max_other_iperp_,
 				       *min_ipar_used_, *max_ipar_used_,
 				       *other_par_,
 				       interp_handle, interp_par_table_handle);
@@ -490,17 +492,12 @@ assert( my_patch() == other_patch()
 //******************************************************************************
 
 //
-// This function "synchronizes" a ghost zone, i.e. it updates the
-// ghost-zone values of the specified gridfns via the appropriate
-// interpatch interpolations.
-//
-// The flags specify which part(s) of the ghost zone we want, but
-// the present implementation only supports the case where all the
-// flags are  true , i.e. we want the entire ghost zone.
+// This function verifies (no-op if ok, error_exit() if not) that the
+// caller wants to operate on the entire ghost zone, since our implementations
+// of synchronize() and compute_Jacobian() only support this case.
 //
-void interpatch_ghost_zone::synchronize
-	(int ghosted_min_gfn, int ghosted_max_gfn,
-	 bool want_min_par_corner,
+void interpatch_ghost_zone::verify_caller_wants_all_of_ghost_zone
+	(bool want_min_par_corner,
 	 bool want_non_corner,
 	 bool want_max_par_corner)
 	const
@@ -508,17 +505,39 @@ void interpatch_ghost_zone::synchronize
 // make sure the caller wants the entire ghost zone
 if (! (want_min_par_corner && want_non_corner && want_max_par_corner))
    then error_exit(ERROR_EXIT,
-"***** interpatch_ghost_zone::synchronize(ghosted_[min,max]_gfn=[%d,%d]):\n"
-"        can only handle full ghost zone, not partial,\n"
-"        but not all flags are true!\n"
+"***** interpatch_ghost_zone::verify_caller_wants_all_of_ghost_zone():\n"
+"        our implementations of  synchronize()  and  compute_Jacobian()\n"
+"        only support operating on the *entire* ghost zone,\n"
+"        but we were passed flags specifying a proper subset!\n"
 "        want_min_par_corner=(int)%d\n"
 "        want_non_corner=(int)%d\n"
 "        want_max_par_corner=(int)%d\n"
 ,
-		   ghosted_min_gfn, ghosted_max_gfn,
 		   want_min_par_corner,
 		   want_non_corner,
 		   want_max_par_corner);			/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
+// This function "synchronizes" a ghost zone, i.e. it updates the
+// ghost-zone values of the specified gridfns via the appropriate
+// interpatch interpolations.
+//
+// The flags specify which part(s) of the ghost zone we want, but
+// the present implementation only supports the case where all the
+// flags are  true , i.e. we want the entire ghost zone.
+//
+void interpatch_ghost_zone::synchronize
+	(int ghosted_min_gfn, int ghosted_max_gfn,
+	 bool want_min_par_corner,
+	 bool want_non_corner,
+	 bool want_max_par_corner)
+{
+verify_caller_wants_all_of_ghost_zone(want_min_par_corner,
+				      want_non_corner,
+				      want_max_par_corner);
 
 // do the interpolation into our result buffer
 other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
@@ -543,3 +562,45 @@ other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
 	}
 	}
 }
+
+//******************************************************************************
+
+//
+// This function allocates the internal buffers for the Jacobian, and
+// computes that Jacobian
+//	    partial synchronize gridfn(ghosted_gfn, iperp, ipar)
+//	------------------------------------------------------------
+//	partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m)
+// where
+//	oiperp = Jacobian_oiperp(iperp)
+//	posn = Jacobian_oipar_posn(iperp, ipar)
+// into the internal buffers.
+//
+void interpatch_ghost_zone::compute_Jacobian
+	(int ghosted_min_gfn, int ghosted_max_gfn,
+	 bool want_min_par_corner,
+	 bool want_non_corner,
+	 bool want_max_par_corner)
+{
+verify_caller_wants_all_of_ghost_zone(want_min_par_corner,
+				      want_non_corner,
+				      want_max_par_corner);
+assert(other_patch_interp_ != NULL);
+
+other_patch_interp_->verify_Jacobian_sparsity_pattern_ok();
+
+other_patch_interp_->molecule_minmax_ipar_m(min_oipar_m_, max_oipar_m_);
+
+if (Jacobian_oipar_posn_ == NULL)
+   then Jacobian_oipar_posn_ = new jtutil::array2d<CCTK_INT>
+				     (min_other_iperp_, max_other_iperp_,
+				      extreme_min_ipar_, extreme_max_ipar_);
+other_patch_interp_->molecule_posn(*Jacobian_oipar_posn_);
+
+if (Jacobian_buffer_ == NULL)
+   then Jacobian_buffer_ = new jtutil::array3d<fp>
+				 (min_other_iperp_, max_other_iperp_,
+				  extreme_min_ipar_, extreme_max_ipar_,
+				  min_oipar_m_, max_oipar_m_);
+other_patch_interp_->Jacobian(*Jacobian_buffer_);
+}