Adding generic types.

src/FiniteDifference/Upwind_gradient.hpp

@@ -35,9 +35,10 @@
  * @param sign: Sign of the velocity with which the wave front is moving.
  * @return Scalar double upwind gradient approximation in the dimension given.
  */
-static double upwinding(double dplus, double dminus, int sign)
+template <typename field_type>
+static field_type upwinding(field_type dplus, field_type dminus, int sign)
 {
-	double grad_upwind = 0;
+	field_type grad_upwind = 0.0;
 	if (dplus * sign < 0
 			&& (dminus + dplus) * sign < 0) grad_upwind = dplus;
 	
@@ -66,9 +67,12 @@ static double upwinding(double dplus, double dminus, int sign)
  * @return Upwind finite difference in one dimension of the property under index Field on the current node with index key.
  */
 template <size_t Field, size_t Sign, typename gridtype, typename keytype>
-double FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
+auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
 {
-	double dplus = 0, dminus = 0;
+//	std::cout << boost::typeindex::type_id_with_cvr<std::decay_t<decltype(first_node)>>() << std::endl;
+	
+	typedef typename std::decay_t<decltype(grid.template get<Field>(key))> field_type;
+	field_type dplus = 0.0, dminus = 0.0;
 	int sign_phi0 = grid.template get<Sign> (key);
 	
 	switch(order)

src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp

@@ -295,4 +295,22 @@ void get_vector_magnitude(gridtype & grid)
 	}
 }
 
+/**@brief Computes the magnitude of the gradient (L2-norm of gradient vector).
+ *
+ * @tparam Phi_grad_in Size_t index of property that contains the gradient.
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the magnitude of gradient should be computed.
+ */
+template <size_t Vector_in, typename magnitude_type, typename key_type, typename gridtype>
+magnitude_type get_vector_magnitude(gridtype & grid, key_type & key)
+{
+	magnitude_type sum = 0;
+	for(size_t d = 0; d < gridtype::dims; d++)
+	{
+		sum += grid.template get<Vector_in> (key)[d] * grid.template get<Vector_in> (key)[d];
+	}
+	return sqrt(sum);
+}
+
+
 #endif //REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP

src/level_set/redistancing_Sussman/RedistancingSussman.hpp

@@ -136,15 +136,16 @@ struct Redist_options
  *
  * @see RedistancingSussman::get_residual_and_change_NB()
  */
+template <typename phi_type>
 struct DistFromSol
 {
-	auto change;   	///< Variable that contains the absolute value of the change of \a &phi; between the
+	phi_type change;   	///< Variable that contains the absolute value of the change of \a &phi; between the
 	///< current
 	///< iteration \a i and the previous iteration \a i-1: @f[ abs( \phi_i - \phi_{i-1} ) @f] It is
 	///< computed for all grid
 	///< points that lie within the narrow band and normalized to the number of grid points that
 	///< lie in that narrow band.
-	auto residual; 	///< Variable that contains the absolute value of how far the gradient magnitude of
+	phi_type residual; 	///< Variable that contains the absolute value of how far the gradient magnitude of
 	///< the current \a &phi; of iteration number \a i is away from being equal to 1: @f[ abs(|\nabla\phi_i| - 1 ) @f]
 	///< It is computed for all grid points that lie within the narrow band and
 	///< normalized to the number of grid points that lie in that narrow band.
@@ -169,7 +170,7 @@ public:
 	 * @param redistOptions User defined options for the Sussman redistancing process
 	 *
 	 */
-	RedistancingSussman(grid_in_type &grid_in, Redist_options &redistOptions) : redistOptions(redistOptions),
+	RedistancingSussman(grid_in_type &grid_in, Redist_options<phi_type> &redistOptions) : redistOptions(redistOptions),
 	                                                                            r_grid_in(grid_in),
 	                                                                            g_temp(grid_in.getDecomposition(),
 	                                                                                   grid_in.getGridInfoVoid().getSize(),
@@ -189,7 +190,7 @@ public:
 	 * gradient of Phi_{n+1},
 	 * sign of the original input Phi_0 (for the upwinding).
 	 */
-	typedef aggregate<phi_type, Point<grid_in_type::dims, phi_type>, int>
+	typedef aggregate<phi_type, phi_type[grid_in_type::dims], int>
 	        props_temp;
 	/** @brief Type definition for the temporary grid.
 	 */
@@ -211,11 +212,12 @@ public:
 	 * having more properties. Computes the gradients. Runs the redistancing on the internal temporary grid. Copies
 	 * resulting signed distance function to the Phi_SDF_out property of the input grid.
 	 */
-	template<size_t Phi_0_in, size_t Phi_SDF_out> void run_redistancing()
+	template <size_t Phi_0_in, size_t Phi_SDF_out>
+	void run_redistancing()
 	{
 		init_temp_grid<Phi_0_in>();
-		init_sign_prop<Phi_n_temp, Phi_0_sign_temp>(
-				g_temp); // initialize Phi_0_sign_temp with the sign of the initial (pre-redistancing) Phi_0
+		init_sign_prop<Phi_n_temp, Phi_0_sign_temp>(g_temp); // Initialize Phi_0_sign_temp with the sign of the
+		// initial (pre-redistancing) Phi_0
 		// Get initial gradients
 		get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, order_upwind_gradient, true);
 		iterative_redistancing(g_temp); // Do the redistancing on the temporary grid
@@ -271,10 +273,11 @@ private:
 	static constexpr size_t Phi_0_sign_temp     = 2; ///< Property index of sign of initial (input) Phi_0 (temp. grid).
 	
 	//	Member variables
-	Redist_options redistOptions; ///< Instantiate redistancing options.
+	Redist_options<phi_type> redistOptions; ///< Instantiate redistancing options.
 	grid_in_type &r_grid_in; ///< Define reference to input grid.
 	
-	DistFromSol distFromSol; ///< Instantiate distance from solution in terms of change, residual, numb. point in NB.
+	DistFromSol<phi_type> distFromSol; ///< Instantiate distance from solution in terms of change, residual, numb.
+	// point in NB.
 	int final_iter = 0; ///< Will be set to the final iteration when redistancing ends.
 	
 	/// Transform the half-bandwidth in no_of_grid_points into physical half-bandwidth kappa.
@@ -282,7 +285,7 @@ private:
 	/**@brief Artificial timestep for the redistancing iterations.
 	 * @see get_time_step_CFL(g_temp_type &grid), get_time_step(), set_user_time_step()
 	 */
-	grid_type::stype time_step;
+	typename grid_in_type::stype time_step;
 	int order_upwind_gradient;
 	//	Member functions
 #ifdef SE_CLASS1
@@ -306,7 +309,7 @@ private:
 	template<size_t Phi_0_in>
 	void init_temp_grid()
 	{
-		phi_type min_value = get_min_val<Phi_0_in>(r_grid_in); // get minimum Phi_0 value on the input grid
+		phi_type min_value = get_min_val<Phi_0_in, phi_type>(r_grid_in); // get minimum Phi_0 value on the input grid
 		init_grid_and_ghost<Phi_n_temp>(g_temp, min_value); // init. Phi_n_temp (incl. ghost) with min. Phi_0
 		copy_gridTogrid<Phi_0_in, Phi_n_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_n_temp
 	}
@@ -321,7 +324,8 @@ private:
 	 * @return Phi_n+1 which is the Phi of the next time step on current node.
 	 *
 	 */
-	phi_type get_phi_nplus1(phi_type phi_n, phi_type phi_n_magnOfGrad, grid_in_type::stype dt, phi_type sgn_phi_n)
+	phi_type get_phi_nplus1(phi_type phi_n, phi_type phi_n_magnOfGrad, typename grid_in_type::stype dt, phi_type
+	sgn_phi_n)
 	{
 		return phi_n + dt * sgn_phi_n * (1 - phi_n_magnOfGrad);
 	}
@@ -339,7 +343,7 @@ private:
 		{
 			auto key = dom.get();
 			const phi_type phi_n = grid.template get<Phi_n_temp>(key);
-			const phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+			const phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp, phi_type>(grid, key);
 			phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
 			grid.template get<Phi_n_temp>(key) = get_phi_nplus1(phi_n, phi_n_magnOfGrad, time_step,
 			                                                         smooth_S(phi_n, epsilon));
@@ -375,7 +379,7 @@ private:
 			if (lays_inside_NB(grid.template get<Phi_n_temp>(key)))
 			{
 				count++;
-				phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+				phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp, phi_type>(grid, key);
 				phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
 				phi_type phi_nplus1 = get_phi_nplus1(grid.template get<Phi_n_temp>(key), phi_n_magnOfGrad, time_step,
 				                                   smooth_S(grid.template get<Phi_n_temp>(key), epsilon));

src/level_set/redistancing_Sussman/tests/convergence_test.cpp

@@ -14,24 +14,27 @@
 #include "l_norms/LNorms.hpp"
 #include "analytical_SDF/AnalyticalSDF.hpp"
 
-double get_min_required_time_step(size_t Nmax, double b_low, double b_up, int dims)
-{
-	double dx = (b_up - b_low) / (double) Nmax;
-	double sum = 0;
-	for (int d = 0; d < dims; d++)
-	{
-		sum += 1 / (dx * dx);
-	}
-	return 0.5 / sum;
-}
+
 
 BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
+	typedef double phi_type;
+	typedef double space_type;
+	const phi_type EPSILON = std::numeric_limits<phi_type>::epsilon();
+	
+	space_type get_min_required_time_step(size_t Nmax, space_type b_low, space_type b_up, int dims)
+	{
+		space_type dx = (b_up - b_low) / (space_type) Nmax;
+		space_type sum = 0;
+		for (int d = 0; d < dims; d++)
+		{
+			sum += 1 / (dx * dx);
+		}
+		return 0.5 / sum;
+	}
 
 	BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_1D)
 	{
 		// CFL dt in 1D for N=64 and c=1 is 0.000503905
-		const double EPSILON = std::numeric_limits<double>::epsilon();
-		std::cout << "epsilon = " << EPSILON << std::endl;
 		
 		const size_t grid_dim = 1;
 		// some indices
@@ -42,17 +45,17 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 		const size_t SDF_exact_grid = 2;
 		const size_t Error_grid = 3;
 		
-		const double box_lower = -1.0;
-		const double box_upper = 1.0;
-		Box<grid_dim, double> box({box_lower}, {box_upper});
+		const space_type box_lower = -1.0;
+		const space_type box_upper = 1.0;
+		Box<grid_dim, space_type> box({box_lower}, {box_upper});
 		Ghost<grid_dim, long int> ghost(0);
-		typedef aggregate<double, double, double, double> props;
-		typedef grid_dist_id<grid_dim, double, props> grid_in_type;
+		typedef aggregate<phi_type, phi_type, phi_type, phi_type> props;
+		typedef grid_dist_id<grid_dim, space_type, props> grid_in_type;
 		
-//		double dt = 0.000503905;
+//		space_type dt = 0.000503905;
 		size_t Nmin = 32;
 		size_t Nmax = 4096;
-		double dt = get_min_required_time_step(Nmax, box_lower, box_upper, grid_dim);
+		space_type dt = get_min_required_time_step(Nmax, box_lower, box_upper, grid_dim);
 		for (size_t N = Nmin; N <= Nmax; N*=2)
 		{
 			const size_t sz[grid_dim] = {N};
@@ -64,14 +67,14 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 			while(dom.isNext())
 			{
 				auto key = dom.get();
-				Point<grid_dim, double> coords = g_dist.getPos(key);
+				Point<grid_dim, space_type> coords = g_dist.getPos(key);
 				g_dist.template get<Phi_0_grid>(key) = sgn(coords.get(x));
 				g_dist.template get<SDF_exact_grid>(key) = coords.get(x);
 				++dom;
 			}
 			
 			{
-				Redist_options redist_options;
+				Redist_options<phi_type> redist_options;
 				redist_options.min_iter = 1e3;
 				redist_options.max_iter = 1e16;
 				
@@ -89,7 +92,7 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				redist_options.print_current_iterChangeResidual = true;     // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF (default: false)
 				redist_options.print_steadyState_iter = true;     // if true, prints out the final iteration number when steady state was reached + final change + residual (default: true)
 				
-				RedistancingSussman<grid_in_type> redist_obj(g_dist,
+				RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist,
 				                                             redist_options);   // Instantiation of Sussman-redistancing class
 				std::cout << "CFL dt = " << redist_obj.get_time_step() << std::endl;
 				redist_obj.set_user_time_step(dt);
@@ -99,14 +102,14 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				
 				// Compute the absolute error between analytical and numerical solution at each grid point
 				get_absolute_error<SDF_sussman_grid, SDF_exact_grid, Error_grid>(g_dist);
-				g_dist.write("grid_postRedistancing", FORMAT_BINARY);
+//				g_dist.write("grid_postRedistancing", FORMAT_BINARY);
 				/////////////////////////////////////////////////////////////////////////////////////////////
 				//	Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
 				//	interface)
 				size_t bc[grid_dim] = {NON_PERIODIC};
-				typedef aggregate<double> props_nb;
-				typedef vector_dist<grid_dim, double, props_nb> vd_type;
-				Ghost<grid_dim, double> ghost_vd(0);
+				typedef aggregate<phi_type> props_nb;
+				typedef vector_dist<grid_dim, space_type, props_nb> vd_type;
+				Ghost<grid_dim, space_type> ghost_vd(0);
 				vd_type vd_narrow_band(0, box, bc, ghost_vd);
 				vd_narrow_band.setPropNames({"error"});
 				const size_t Error_vd = 0;
@@ -115,7 +118,7 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				NarrowBand<grid_in_type> narrowBand_points(g_dist, narrow_band_width); // Instantiation of NarrowBand class
 				narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
 				                                                                                                 vd_narrow_band);
-			    vd_narrow_band.write("vd_error_N" + std::to_string(N), FORMAT_BINARY);
+//			    vd_narrow_band.write("vd_error_N" + std::to_string(N), FORMAT_BINARY);
 //				vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
 				// Compute the L_2- and L_infinity-norm and save to file
 				L_norms lNorms_vd;
@@ -128,8 +131,6 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 	
 	BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_3D)
 	{
-		const double EPSILON = std::numeric_limits<double>::epsilon();
-	
 		const size_t grid_dim = 3;
 		// some indices
 		const size_t x                      = 0;
@@ -143,24 +144,24 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 		
 		for (size_t N=32; N <=128; N*=2)
 		{
-			const double dt = 0.000165334;
+			const space_type dt = 0.000165334;
 			const size_t sz[grid_dim] = {N, N, N};
-			const double radius = 1.0;
-			const double box_lower = 0.0;
-			const double box_upper = 4.0 * radius;
-			Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
+			const space_type radius = 1.0;
+			const space_type box_lower = 0.0;
+			const space_type box_upper = 4.0 * radius;
+			Box<grid_dim, space_type> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
 			Ghost<grid_dim, long int> ghost(0);
-			typedef aggregate<double, double, double, double> props;
-			typedef grid_dist_id<grid_dim, double, props > grid_in_type;
+			typedef aggregate<phi_type, phi_type, phi_type, phi_type> props;
+			typedef grid_dist_id<grid_dim, space_type, props > grid_in_type;
 			grid_in_type g_dist(sz, box, ghost);
 			g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
 			
-			const double center[grid_dim] = {0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower)};
+			const space_type center[grid_dim] = {0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower)};
 			init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
 			
 			int order = 1;
 			{
-				Redist_options redist_options;
+				Redist_options<phi_type> redist_options;
 				redist_options.min_iter                             = 1e4;
 				redist_options.max_iter                             = 1e4;
 				
@@ -174,7 +175,8 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				redist_options.print_current_iterChangeResidual     = true;     // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF (default: false)
 				redist_options.print_steadyState_iter               = true;     // if true, prints out the final iteration number when steady state was reached + final change + residual (default: true)
 				
-				RedistancingSussman<grid_in_type> redist_obj(g_dist, redist_options);   // Instantiation of Sussman-redistancing class
+				RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options);   // Instantiation of
+				// Sussman-redistancing class
 				redist_obj.set_user_time_step(dt);
 				std::cout << "dt set to = " << dt << std::endl;
 				// Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
@@ -191,9 +193,9 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				//	Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
 				//	interface)
 				size_t bc[grid_dim] = {PERIODIC, PERIODIC, PERIODIC};
-				typedef aggregate<double> props_nb;
-				typedef vector_dist<grid_dim, double, props_nb> vd_type;
-				Ghost<grid_dim, double> ghost_vd(0);
+				typedef aggregate<phi_type> props_nb;
+				typedef vector_dist<grid_dim, space_type, props_nb> vd_type;
+				Ghost<grid_dim, space_type> ghost_vd(0);
 				vd_type vd_narrow_band(0, box, bc, ghost_vd);
 				vd_narrow_band.setPropNames({"error"});
 				const size_t Error_vd = 0;
@@ -202,8 +204,8 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				NarrowBand<grid_in_type> narrowBand_points(g_dist, narrow_band_width); // Instantiation of NarrowBand class
 				narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
 				                                                                                                 vd_narrow_band);
-				vd_narrow_band.write("vd_nb8p_error_N" + std::to_string(N) + "_order" +
-				std::to_string(order), FORMAT_BINARY);
+//				vd_narrow_band.write("vd_nb8p_error_N" + std::to_string(N) + "_order" +
+//				std::to_string(order), FORMAT_BINARY);
 //				vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
 				// Compute the L_2- and L_infinity-norm and save to file
 				L_norms lNorms_vd;

src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp

@@ -51,7 +51,7 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
 		init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
 		
 		
-		Redist_options redist_options;
+		Redist_options<phi_type> redist_options;
 		redist_options.min_iter                             = 1e4;
 		redist_options.max_iter                             = 1e4;