Commit 0e987145 authored by Mario Morales Hernandez's avatar Mario Morales Hernandez
Browse files

curvilinear bc implementation. not tested yet. 

This update introduces support for curvilinear boundary conditions,
allowing the simulation to handle complex geometries that do not
align with the Cartesian axes. The implementation adds functionality
to the ExtBC class to load boundary polylines and compute local normal
vectors, which are then used in the compute_extbc_values kernel to
project discharge components correctly. To maintain consistency with
the underlying grid topology, a dominant-component logic was
implemented to ensure ghost cell indexing always utilizes
face-connected neighbors rather than invalid diagonal offsets.
Additionally, the matrix class was extended with a new member function
to handle DEM ghost cell updates for curvilinear boundaries, ensuring
proper data encapsulation and robust initialization across the domain.
parent ec35bc4b

src/config_utils.h

+110 −19
Original line number Diff line number Diff line
@@ -91,7 +91,8 @@ namespace ConfigUtils 


		std::vector<int> extbc_bctype;	/**< Contains all external boundary condition type serially. */

		std::vector<std::string> extbc_fname;	/**< Contains all external boundary condition file name serially. */



		std::vector<std::string> extbc_polyline_fname;  // Polyline file for each BC

		std::vector<bool> extbc_is_curvilinear;         // Flag for curvilinear BCs

	};
@@ -133,6 +134,17 @@ namespace ConfigUtils 


/** @brief It extracts each external boundary condition file and constructs an attribute key-value mapping.

*

*  External BC file format:

*  - Axis-aligned: geometry_flag, bc_type, x1, y1, x2, y2[, bc_data_file]

*    geometry_flag = 0 for axis-aligned boundaries

*    bc_type: 0=open boundaries (zero gradient), 1=level time-varying , 2=normal slope, 3=froude number

*    bc_data_file: required for bc_type=1, otherwise a number representing fixed value

*  - Curvilinear: geometry_flag, bc_type, polyline_file, bc_data_file

*    geometry_flag = 1 for curvilinear boundaries

*    bc_type: 0=open boundaries (zero gradient), 1=level time-varying , 2=normal slope, 3=froude number

*    polyline_file: file containing boundary geometry

*    bc_data_file: required for bc_type=1, otherwise a number representing fixed value

*

*  @param filename file to parse

*  @param dir parent directory of filename

*  @return Attribute key value mapping
@@ -321,6 +333,8 @@ namespace ConfigUtils 
		std::string src_x1_loc = "", src_y1_loc = "";

		std::string src_x2_loc = "", src_y2_loc = "";

		std::string bcfname = "";

		std::string polyline_fname = "";

		std::string is_curvilinear = "";



		while (std::getline(ifs, line))

		{
@@ -334,22 +348,64 @@ namespace ConfigUtils 
					line.erase(line.begin() + incomm);

				}

				std::vector<std::string> kv = StringUtils::split(line, ',');

				std::string bctype_value = kv[0];

				std::string x1_value = kv[1];

				std::string y1_value = kv[2];

				std::string x2_value = kv[3];

				std::string y2_value = kv[4];

				

				int auxtype=std::stoi(kv[0]);

				if (kv.size() < 3) {

					std::cerr << ERROR "Invalid external BC format. Expected at least 3 fields." << std::endl;

					exit(EXIT_FAILURE);

				}

				

				int geometry_flag = std::stoi(kv[0]);

				std::string bctype_value = kv[1];

				int auxtype = std::stoi(kv[1]);

				

				std::string x1_value, y1_value, x2_value, y2_value;

				std::string bcfname_value = "";

				std::string polyline_value = "";

				std::string curvilinear_flag = (geometry_flag == 1) ? "1" : "0";

				

				if (geometry_flag == 1) {

					if (kv.size() < 4) {

						std::cerr << ERROR "Curvilinear BC requires: geometry_flag, bc_type, polyline_file, bc_data_file" << std::endl;

						exit(EXIT_FAILURE);

					}

					polyline_value = kv[2];

					bcfname_value = kv[3];

					x1_value = "0";

					y1_value = "0";

					x2_value = "0";

					y2_value = "0";

				} else if (geometry_flag == 0) {

					if (kv.size() < 6) {

						std::cerr << ERROR "Axis-aligned BC requires: geometry_flag, bc_type, x1, y1, x2, y2" << std::endl;

						exit(EXIT_FAILURE);

					}

					x1_value = kv[2];

					y1_value = kv[3];

					x2_value = kv[4];

					y2_value = kv[5];

					polyline_value = "";

					

					if(auxtype != 0)

					{

					bcfname_value = kv[5];

						if (kv.size() >= 7) {

							bcfname_value = kv[6];

						} else {

							std::cerr << ERROR "BC type " << auxtype << " requires a data file" << std::endl;

							exit(EXIT_FAILURE);

						}

					}

					else

					{

						if (kv.size() >= 7) {

							bcfname_value = kv[6];

						} else {

							bcfname_value = "0.0";

						}

					}

				} else {

					std::cerr << ERROR "Invalid geometry flag: " << geometry_flag << ". Use 0 for axis-aligned, 1 for curvilinear." << std::endl;

					exit(EXIT_FAILURE);

				}



				if (bctype_value[0] == '"' && bctype_value[bctype_value.size() - 1] == '"')

				{
@@ -385,22 +441,45 @@ namespace ConfigUtils 
					bcfname_value.erase(bcfname_value.end() - 1);

				}

				

				if (!polyline_value.empty() && polyline_value[0] == '"' && polyline_value[polyline_value.size() - 1] == '"')

				{

					polyline_value.erase(polyline_value.begin());

					polyline_value.erase(polyline_value.end() - 1);

				}



				bctype = bctype + "," + bctype_value;

				src_x1_loc = src_x1_loc + "," + x1_value;

				src_y1_loc = src_y1_loc + "," + y1_value;

				src_x2_loc = src_x2_loc + "," + x2_value;

				src_y2_loc = src_y2_loc + "," + y2_value;

				is_curvilinear = is_curvilinear + "," + curvilinear_flag;



				if(geometry_flag == 1)

				{

					polyline_fname = polyline_fname + "," + dir + "/" + polyline_value;

					if(auxtype == 1)

					{

						bcfname = bcfname + "," + dir + "/" + bcfname_value;

					}

					else

					{

						bcfname = bcfname + "," + bcfname_value;

					}

				}

				else

				{

					polyline_fname = polyline_fname + ",";

					if(auxtype == 1)

					{

						bcfname = bcfname + "," + dir + "/" + bcfname_value;

					}

				else   //case 0, 2, 3

					else

					{

						bcfname = bcfname + "," + bcfname_value;

					}

				}

			}

		}



		bctype = bctype.substr(1);

		src_x1_loc = src_x1_loc.substr(1);
@@ -408,6 +487,8 @@ namespace ConfigUtils 
		src_x2_loc = src_x2_loc.substr(1);

		src_y2_loc = src_y2_loc.substr(1);

		bcfname = bcfname.substr(1);

		polyline_fname = polyline_fname.substr(1);

		is_curvilinear = is_curvilinear.substr(1);



		arglist.insert(std::pair<std::string, std::string>("extbc_bctype", bctype));

		arglist.insert(std::pair<std::string, std::string>("extbc_x1_loc", src_x1_loc));
@@ -415,9 +496,12 @@ namespace ConfigUtils 
		arglist.insert(std::pair<std::string, std::string>("extbc_x2_loc", src_x2_loc));

		arglist.insert(std::pair<std::string, std::string>("extbc_y2_loc", src_y2_loc));

		arglist.insert(std::pair<std::string, std::string>("extbc_fname", bcfname));

		arglist.insert(std::pair<std::string, std::string>("extbc_polyline_fname", polyline_fname));

		arglist.insert(std::pair<std::string, std::string>("extbc_is_curvilinear", is_curvilinear));



		ifs.close();





		return arglist;

	}
@@ -506,6 +590,13 @@ namespace ConfigUtils 
			arglist.extbc_x2_loc = StringUtils::vecstr_to_vecflt<T>(StringUtils::split((args("extbc_x2_loc", extbc_map)), ','));

			arglist.extbc_y2_loc = StringUtils::vecstr_to_vecflt<T>(StringUtils::split((args("extbc_y2_loc", extbc_map)), ','));

			arglist.extbc_fname = StringUtils::split((args("extbc_fname", extbc_map)), ',');

			arglist.extbc_polyline_fname = StringUtils::split((args("extbc_polyline_fname", extbc_map)), ',');



			std::vector<int> curv_flags = StringUtils::vecstr_to_vecint(StringUtils::split((args("extbc_is_curvilinear", extbc_map)), ','));

			arglist.extbc_is_curvilinear.resize(curv_flags.size());

			for (size_t i = 0; i < curv_flags.size(); i++) {

				arglist.extbc_is_curvilinear[i] = (curv_flags[i] != 0);

			}

		}

		

		if(strcmp(arglist.input_option.c_str(), "")==0){ //no input_option provided (for legacy version)

src/constants.h

+14 −10
Original line number Diff line number Diff line
@@ -102,6 +102,8 @@ typedef double value_t; /**< Data type to represent floating-point number. It 
#define OBSH 20    /**< Water depth for observation point array position in vector. */

#define OBSQX 21    /**< Flux X for observation point array position in vector. */

#define OBSQY 22    /**< Flux Y for observation point array position in vector. */

#define BCNORMALX 23    /**< Boundary condition normal X component array position in vector. */

#define BCNORMALY 24    /**< Boundary condition normal Y component array position in vector. */



#define SRCP 0    /**< Flow locations index array position in vector. */

#define RUNID 1    /**< Runoff id array position in vector. */
@@ -110,6 +112,8 @@ typedef double value_t; /**< Data type to represent floating-point number. It 
#define BCINDEXSTART 4    /**< Boundary condition's start index array position in vector. */

#define BCNROWSVARS 5    /**< Boundary condition's number of rows variable array position in vector. */

#define OBSRELATIVEINDEX 6    /**< Observation point index array after domain decomposition position in vector. */

#define BCCOLS 7    /**< Boundary condition cell columns array position in vector. */

#define BCROWS 8    /**< Boundary condition cell rows array position in vector. */



#define TIMER_NSECS 0    /**< To use nano second in Timer. */

#define TIMER_SECS 1    /**< To use second in Timer. */

src/extbc.h

+599 −341
Original line number Diff line number Diff line 
/** @file extbc.h

 *  @brief Header containing the ExtBC class

 *  @brief Header containing the unified ExtBC class for both axis-aligned and curvilinear boundaries

 *

 *  This contains the subroutines and eventually any 

 *  macros, constants, etc. needed for ExtBC class

 *  Supports both axis-aligned and curvilinear (polyline-based) external boundary conditions

 *

 *  @author Mario Morales Hernandez

 *  @author Sudershan Gangrade
@@ -19,6 +20,13 @@ 
#ifndef EXTBC_H

#define EXTBC_H



#include <vector>

#include <string>

#include <cmath>

#include <algorithm>

#include <iostream>

#include <fstream>

#include <sstream>



namespace ExtBC

{
@@ -32,13 +40,18 @@ namespace ExtBC 
*/			

		extBC();



		

/** @brief Constructor. Takes filename containing boundary condition and boundary condition type. Reads from files and push each row in a vector and construct data.

*

*  @param filename File name

*  @param bctype Boundary condition type

		/**

		 * @brief Constructor for external boundary condition

		 * @param filename_or_data For bctype=1 (time-varying): filename of BC data file.

		 *                         For bctype=2,3 (constant): string representation of constant value

		 * @param bctype Hydraulic boundary condition type:

		 *               0 = open boundaries (zero gradient)

		 *               1 = time-varying level (from file)

		 *               2 = normal slope

		 *               3 = froude number

		 * @param is_curvilinear Geometry flag: true for curvilinear (polyline-based), false for axis-aligned

		 */

		extBC(std::string filename, int bctype);

		extBC(std::string filename_or_data, int bctype, bool is_curvilinear = false);

		



/** @brief Reads from files and push each row in a vector and construct data.
@@ -48,6 +61,7 @@ namespace ExtBC 
*/			

void load_from_file(std::string filename, int bctype);



		void load_polyline_from_file(std::string polyline_filename);

				

/** @brief It checks for extreme boundary condition and calculates the number of cells in that boundary condition.

*
@@ -70,11 +84,16 @@ namespace ExtBC 
*/	

		void create_involved_cells(std::vector<int> e_cols, std::vector<int> e_rows, int ncols, int nrows, int bctype);

		

		int compute_cells_along_polyline(const std::vector<T>& polyline_x, 

		                                  const std::vector<T>& polyline_y,

		                                  int ncols, int nrows,

		                                  T dx, T dy, T xllcorner, T yllcorner);

		

		void create_involved_cells_curvilinear(const std::vector<T>& polyline_x,

		                                        const std::vector<T>& polyline_y,

		                                        int ncols, int nrows,

		                                        T dx, T dy, T xllcorner, T yllcorner);

		

/** @brief It returns all data saved from boundary condition file.

*

*  @return Boundary condition data

*/

		std::vector<std::vector<T>> get_rows();
@@ -112,18 +131,31 @@ namespace ExtBC 
*/

		void convert_to_secs();

		

		

		int ncells;	/**< Number of cells of a boundary condition. */

		int location;	/**< 0--> westBoundary  1-->northBoundary 2-->eastBoundary 3--> southBoundary */

		int ncells_local;	/**< Number of cells of a boundary condition in a subdomain. */

		std::vector<int> extreme_rows;	/**< Extreme rows */

		std::vector<int>extreme_cols;	/**< Extreme columns */

		std::vector<int> i_cols;	/**< Involved columns. */

		std::vector<int> i_rows;	/**< Involved rows. */

		void compute_boundary_normals(const std::vector<T>& polyline_x,

		                               const std::vector<T>& polyline_y,

		                               const std::vector<int>& i_cols,

		                               const std::vector<int>& i_rows,

		                               std::vector<T>& normal_x,

		                               std::vector<T>& normal_y,

		                               T dx, T dy, T xllcorner, T yllcorner,

		                               int ncols, int nrows);

		

		bool is_curvilinear_;

		int ncells;

		int location;

		int ncells_local;

		std::vector<int> extreme_rows;

		std::vector<int> extreme_cols;

		std::vector<int> i_cols;

		std::vector<int> i_rows;

		std::vector<T> polyline_x_;

		std::vector<T> polyline_y_;

		std::vector<T> normal_x_;

		std::vector<T> normal_y_;

		

	private:

		int num_rows_;	/**< Number of rows. */

		std::vector<std::vector<T> > data_;	/**< Data contains all the rows of a boundary condition file. */

		int num_rows_;

		std::vector<std::vector<T> > data_;

	};
@@ -131,16 +163,25 @@ namespace ExtBC 
	extBC<T>::extBC()

	{

		set_num_rows(0);

		is_curvilinear_ = false;

		ncells = 0;

		ncells_local = 0;

		location = -1;

	}





	template<class T>

	extBC<T>::extBC(std::string filename, int bctype)

	extBC<T>::extBC(std::string filename_or_data, int bctype, bool is_curvilinear)

	{

		is_curvilinear_ = is_curvilinear;

		ncells = 0;

		ncells_local = 0;

		location = -1;

		

		if(bctype==1){

			load_from_file(filename, bctype);

			load_from_file(filename_or_data, bctype);

		}else{	

			std::istringstream iss(filename);

			std::istringstream iss(filename_or_data);

			T temp;

			iss >> temp;

			std::vector<T> row;
@@ -201,6 +242,40 @@ namespace ExtBC 
	}





	template<typename T>

	void extBC<T>::load_polyline_from_file(std::string polyline_filename)

	{

		std::ifstream infile(polyline_filename.c_str());

		if (!infile.is_open())

		{

			std::cerr << "ERROR: Cannot open polyline file: " << polyline_filename << std::endl;

			exit(EXIT_FAILURE);

		}



		polyline_x_.clear();

		polyline_y_.clear();



		std::string line;

		while (std::getline(infile, line))

		{

			std::istringstream iss(line);

			T x, y;

			if (iss >> x >> y)

			{

				polyline_x_.push_back(x);

				polyline_y_.push_back(y);

			}

		}

		infile.close();



		if (polyline_x_.size() < 2)

		{

			std::cerr << "ERROR: Polyline must have at least 2 points" << std::endl;

			exit(EXIT_FAILURE);

		}

	}





	template<typename T>

	int extBC<T>::check_extreme_extbc(std::vector<int> e_cols, std::vector<int> e_rows, int ncols, int nrows)

	{
@@ -292,6 +367,100 @@ namespace ExtBC 
	}





	template<typename T>

	int extBC<T>::compute_cells_along_polyline(const std::vector<T>& polyline_x, 

	                                            const std::vector<T>& polyline_y,

	                                            int ncols, int nrows,

	                                            T dx, T dy, T xllcorner, T yllcorner)

	{

		std::vector<std::pair<int, int>> cell_set;

		

		for (size_t seg = 0; seg < polyline_x.size() - 1; seg++)

		{

			T x1 = polyline_x[seg];

			T y1 = polyline_y[seg];

			T x2 = polyline_x[seg + 1];

			T y2 = polyline_y[seg + 1];

			

			T seg_length = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));

			T cell_size = std::min(dx, dy);

			int num_samples = std::max(2, static_cast<int>(seg_length / cell_size * 2.0));

			

			for (int i = 0; i <= num_samples; i++)

			{

				T t = static_cast<T>(i) / static_cast<T>(num_samples);

				T x = x1 + t * (x2 - x1);

				T y = y1 + t * (y2 - y1);

				

				int col = static_cast<int>((x - xllcorner) / dx);

				int row = static_cast<int>((yllcorner + nrows * dy - y) / dy);

				

				if (col >= 0 && col < ncols && row >= 0 && row < nrows)

				{

					std::pair<int, int> cell(col, row);

					if (std::find(cell_set.begin(), cell_set.end(), cell) == cell_set.end())

					{

						cell_set.push_back(cell);

					}

				}

			}

		}

		

		return cell_set.size();

	}





	template<typename T>

	void extBC<T>::create_involved_cells_curvilinear(const std::vector<T>& polyline_x,

	                                                   const std::vector<T>& polyline_y,

	                                                   int ncols, int nrows,

	                                                   T dx, T dy, T xllcorner, T yllcorner)

	{

		std::vector<std::pair<int, int>> cell_set;

		

		for (size_t seg = 0; seg < polyline_x.size() - 1; seg++)

		{

			T x1 = polyline_x[seg];

			T y1 = polyline_y[seg];

			T x2 = polyline_x[seg + 1];

			T y2 = polyline_y[seg + 1];

			

			T seg_length = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));

			T cell_size = std::min(dx, dy);

			int num_samples = std::max(2, static_cast<int>(seg_length / cell_size * 2.0));

			

			for (int i = 0; i <= num_samples; i++)

			{

				T t = static_cast<T>(i) / static_cast<T>(num_samples);

				T x = x1 + t * (x2 - x1);

				T y = y1 + t * (y2 - y1);

				

				int col = static_cast<int>((x - xllcorner) / dx);

				int row = static_cast<int>((yllcorner + nrows * dy - y) / dy);

				

				if (col >= 0 && col < ncols && row >= 0 && row < nrows)

				{

					std::pair<int, int> cell(col, row);

					if (std::find(cell_set.begin(), cell_set.end(), cell) == cell_set.end())

					{

						cell_set.push_back(cell);

					}

				}

			}

		}

		

		ncells = cell_set.size();

		i_cols.resize(ncells);

		i_rows.resize(ncells);

		

		for (int i = 0; i < ncells; i++)

		{

			i_cols[i] = cell_set[i].first;

			i_rows[i] = cell_set[i].second;

		}

	}





	template<typename T>

	std::vector<std::vector<T>> extBC<T>::get_rows()

	{
@@ -336,6 +505,95 @@ namespace ExtBC 
			(*it)[0] = it->at(0) * HOUR_TO_SEC_FACTOR;

		}

	}



	template<typename T>

	void extBC<T>::compute_boundary_normals(const std::vector<T>& polyline_x,

	                               const std::vector<T>& polyline_y,

	                               const std::vector<int>& i_cols,

	                               const std::vector<int>& i_rows,

	                               std::vector<T>& normal_x,

	                               std::vector<T>& normal_y,

	                               T dx, T dy, T xllcorner, T yllcorner,

	                               int ncols, int nrows)

	{

		int ncells = i_cols.size();

		normal_x.resize(ncells);

		normal_y.resize(ncells);

		

		for (int i = 0; i < ncells; i++)

		{

			T cell_x = xllcorner + i_cols[i] * dx + dx * 0.5;

			T cell_y = yllcorner + (nrows - i_rows[i]) * dy - dy * 0.5;

			

			T min_dist = 1e30;

			int closest_seg = 0;

			T closest_t = 0.0;

			

			for (size_t seg = 0; seg < polyline_x.size() - 1; seg++)

			{

				T x1 = polyline_x[seg];

				T y1 = polyline_y[seg];

				T x2 = polyline_x[seg + 1];

				T y2 = polyline_y[seg + 1];

				

				T dx_seg = x2 - x1;

				T dy_seg = y2 - y1;

				T seg_len_sq = dx_seg * dx_seg + dy_seg * dy_seg;

				

				T t = 0.0;

				if (seg_len_sq > 1e-10)

				{

					t = ((cell_x - x1) * dx_seg + (cell_y - y1) * dy_seg) / seg_len_sq;

					t = std::max(T(0.0), std::min(T(1.0), t));

				}

				

				T proj_x = x1 + t * dx_seg;

				T proj_y = y1 + t * dy_seg;

				

				T dist = sqrt((cell_x - proj_x) * (cell_x - proj_x) + 

				              (cell_y - proj_y) * (cell_y - proj_y));

				

				if (dist < min_dist)

				{

					min_dist = dist;

					closest_seg = seg;

					closest_t = t;

				}

			}

			

			T x1 = polyline_x[closest_seg];

			T y1 = polyline_y[closest_seg];

			T x2 = polyline_x[closest_seg + 1];

			T y2 = polyline_y[closest_seg + 1];

			

			T tangent_x = x2 - x1;

			T tangent_y = y2 - y1;

			T tangent_len = sqrt(tangent_x * tangent_x + tangent_y * tangent_y);

			

			if (tangent_len > 1e-10)

			{

				tangent_x /= tangent_len;

				tangent_y /= tangent_len;

			}

			

			T norm_x = -tangent_y;

			T norm_y = tangent_x;

			

			T proj_x = x1 + closest_t * (x2 - x1);

			T proj_y = y1 + closest_t * (y2 - y1);

			T to_cell_x = cell_x - proj_x;

			T to_cell_y = cell_y - proj_y;

			

			if (norm_x * to_cell_x + norm_y * to_cell_y < 0)

			{

				norm_x = -norm_x;

				norm_y = -norm_y;

			}

			

			normal_x[i] = norm_x;

			normal_y[i] = norm_y;

		}

	}

}



#endif

src/kernels.h

+129 −58

File changed.

Preview size limit exceeded, changes collapsed.

src/matrix.h

+64 −1
Original line number Diff line number Diff line
@@ -284,6 +284,21 @@ namespace Matrix 
*/		

		void copy_value_into_ghost_cells_location(std::vector<int> irows, std::vector<int> icols, int ncells, int location);



/** @brief It copies the elevation of boundary cells values into ghost cells for curvilinear grids.

*

*  @param irows Index of boundary cells row

*  @param icols Index of boundary cells column

*  @param normal_x Normal vector x component

*  @param normal_y Normal vector y component

*  @param ncells Number of cells

*/	

		void copy_value_into_ghost_cells_curvilinear(

			std::vector<int> irows, 

			std::vector<int> icols, 

			std::vector<T> normal_x,

			std::vector<T> normal_y,

			int ncells);



		

/** @brief Put infinite walls in boundary cells.

*
@@ -899,6 +914,54 @@ namespace Matrix 
		}

	}



	template<typename T>

	void matrix<T>::copy_value_into_ghost_cells_curvilinear(

		std::vector<int> irows, 

		std::vector<int> icols, 

		std::vector<T> normal_x,

		std::vector<T> normal_y,

		int ncells)

	{

		int ncols = this->cols_;

		

		for (int i = 0; i < ncells; i++)

		{

			int ix = irows[i] + GHOST_CELL_PADDING;

			int iy = icols[i] + GHOST_CELL_PADDING;

			T nx = normal_x[i];

			T ny = normal_y[i];

			

			int ii = ix * (long long)ncols + iy;

			int ib;

			

			// Choose ghost cell based on dominant normal component

			T abs_nx = fabs(nx);

			T abs_ny = fabs(ny);

			

			if (abs_nx >= abs_ny) {

			// Normal predominantly horizontal → ghost is east or west

			if (nx > 0.0) {

				ib = ii + 1;              // east neighbor (column +1)

			} else if (nx < 0.0) {

				ib = ii - 1;              // west neighbor (column -1)

			} else {

				continue;                 // degenerate normal, skip

			}

			} else {

			// Normal predominantly vertical → ghost is north or south

			if (ny > 0.0) {

				ib = ii + ncols;          // south neighbor (row +1)

			} else if (ny < 0.0) {

				ib = ii - ncols;          // north neighbor (row -1)

			} else {

				continue;                 // degenerate normal, skip

			}

			}

			

			// Copy DEM value from interior to ghost cell

			this->data_[ib] = this->data_[ii];

		}

	}



	template<typename T>

	void matrix<T>::set_infinite_walls()