add 1D GridSubdomain and use in 2D and 3D implementations

Author: soypat

md1/md1.go

@@ -162,3 +162,31 @@ func (nra NewtonRaphsonSolver) Root(x0 float64, f func(xGuess float64) float64)
 	}
 	return x_root, -nra.MaxIterations
 }
+
+// GridSubdomain returns the indexing for a 1D subdomain contained within a larger gridded domain of equally spaced nGridPoints points.
+//   - iStart contains the starting index of the first point contained within the subdomain.
+//   - nSub is the number of consecutive grid points contained within the subdomain.
+//
+// This function is primarily used for higher dimension GridSubdomain calls in 2D and 3D packages
+func GridSubdomain(domMin, domMax float64, nGridPoints int, subMin, subMax float64) (iStart, nSub int) {
+	if domMin > domMax || subMin > subMax || subMax > domMax || nGridPoints <= 1 {
+		panic("invalid arguments to gridSubdomain")
+	}
+	dx := (domMax - domMin) / float64(nGridPoints-1) // The size of grid spaces.
+	startOffset := subMin - domMin
+	iStartf := startOffset / dx
+	iStart = min(int(math.Ceil(iStartf)), nGridPoints-1)
+
+	endOffset := subMax - domMin
+	iEndf := endOffset / dx
+	iEnd := min(int(math.Ceil(iEndf)), nGridPoints-1)
+	nSub = iEnd - iStart
+	return iStart, nSub
+}
+
+func min(a, b int) int {
+	if a < b {
+		return a
+	}
+	return b
+}

md1/md1_test.go

@@ -0,0 +1,60 @@
+// DO NOT EDIT.
+// This file was generated automatically
+// from gen.go. Please do not edit this file.
+
+package md1
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func Test_gridSubdomainInternal(t *testing.T) {
+	rng := rand.New(rand.NewSource(1))
+	for i := 0; i < 100000; i++ {
+		dmin := float64(rng.Float64())
+		dmax := dmin + float64(rng.Float64())
+		dsize := dmax - dmin
+
+		smin := dmin + float64(rng.Float64())*dsize
+		smaxsize := dmax - smin
+		smax := smin + float64(rng.Float64())*smaxsize
+		if i%8 == 0 {
+			smax = dmax
+		}
+		if i%16 == 0 {
+			smin = dmin
+		}
+		ssize := smax - smin
+		for Nd := 2; Nd < 10; Nd++ {
+			istart, nsub := GridSubdomain(dmin, dmax, Nd, smin, smax)
+			if nsub > Nd {
+				t.Error("subdomain contains more grid points than actual domain", nsub, Nd)
+			} else if istart >= Nd {
+				t.Error("istart >= Nd", istart, Nd)
+			} else if istart+nsub > Nd {
+				t.Error("istart+nsub>Nd", istart, nsub, istart+nsub, Nd)
+			} else if nsub < 0 {
+				t.Error("nsub<0", nsub)
+			}
+			dx := dsize / float64(Nd-1)
+			calc_xmax := dmin + float64(istart+nsub-1)*dx
+			calc_xmin := dmin + float64(istart)*dx
+			if nsub > 0 && calc_xmax > smax {
+				t.Error("xmax > smax", calc_xmax, smax)
+			}
+			if nsub > 0 && calc_xmin < smin {
+				t.Error("xmin < smin", calc_xmin, smin)
+			}
+			if nsub > 0 && calc_xmax > dmax {
+				t.Error("xmax > dmax", calc_xmax, dmax)
+			}
+			_ = ssize
+			if t.Failed() {
+				istart, nsub = GridSubdomain(dmin, dmax, Nd, smin, smax) // For debugging.
+				t.Fatalf("failed on output: i=%d ns=%d  input: d=(%.3f,%.3f) N=%d sub=(%.3f, %.3f) ix=%.3f",
+					istart, nsub, dmin, dmax, Nd, smin, smax, calc_xmin)
+			}
+		}
+	}
+}

md2/grid.go

@@ -4,9 +4,7 @@
 
 package md2
 
-import (
-	math "math"
-)
+import ms1 "github.com/soypat/geometry/md1"
 
 // AppendGrid splits the argument bounds [Box] x,y axes by nx,ny, respectively
 // and generates points on the vertices generated by the division and appends them to dst, returning the result.
@@ -50,23 +48,8 @@ func GridSubdomain(domain Box, nxDomain, nyDomain int, subdomain Box) (iStart, n
 	if !domain.ContainsBox(subdomain) {
 		panic("subdomain not contained in domain")
 	}
-	dx := (domain.Max.X - domain.Min.X) / float64(nxDomain-1)
-	dy := (domain.Max.Y - domain.Min.Y) / float64(nyDomain-1)
-
-	off := Sub(subdomain.Min, domain.Min)
-	ix0 := iceil(off.X / dx)
-	iy0 := iceil(off.Y / dy)
-	iStart = ix0 + iy0*nxDomain
-
-	offEnd := Sub(subdomain.Max, domain.Min)
-	ixf := int(offEnd.X / dx)
-	iyf := int(offEnd.Y / dy)
-
-	nxSub = ixf - ix0 + 1
-	nySub = iyf - iy0 + 1
-	return iStart, nxSub, nySub
-}
-
-func iceil(f float64) int {
-	return int(math.Ceil(f))
+	var ix0, iy0 int
+	ix0, nxSub = ms1.GridSubdomain(domain.Min.X, domain.Max.X, nxDomain, subdomain.Min.X, subdomain.Max.X)
+	iy0, nySub = ms1.GridSubdomain(domain.Min.Y, domain.Max.Y, nyDomain, subdomain.Min.Y, subdomain.Max.Y)
+	return ix0 + iy0*nxDomain, nxSub, nySub
 }

md2/md2_test.go

@@ -4,7 +4,9 @@
 
 package md3
 
-import math "math"
+import (
+	ms1 "github.com/soypat/geometry/md1"
+)
 
 // AppendGrid splits the argument bounds [Box] x,y,z axes by nx,ny,nz, respectively
 // and generates points on the vertices generated by the division and appends them to dst, returning the result.
@@ -54,26 +56,9 @@ func GridSubdomain(domain Box, nxDomain, nyDomain, nzDomain int, subdomain Box)
 	if !domain.ContainsBox(subdomain) {
 		panic("subdomain not contained in domain")
 	}
-	dx := (domain.Max.X - domain.Min.X) / float64(nxDomain-1)
-	dy := (domain.Max.Y - domain.Min.Y) / float64(nyDomain-1)
-	dz := (domain.Max.Z - domain.Min.Z) / float64(nzDomain-1)
-	off := Sub(subdomain.Min, domain.Min)
-	ix0 := iceil(off.X / dx)
-	iy0 := iceil(off.Y / dy)
-	iz0 := iceil(off.Z / dz)
-	iStart = ix0 + iy0*nxDomain + iz0*(nxDomain+nyDomain)
-
-	offEnd := Sub(subdomain.Max, domain.Min)
-	ixf := int(offEnd.X / dx)
-	iyf := int(offEnd.Y / dy)
-	izf := int(offEnd.Z / dz)
-
-	nxSub = ixf - ix0 + 1
-	nySub = iyf - iy0 + 1
-	nzSub = izf - iz0 + 1
-	return iStart, nxSub, nySub, nzSub
-}
-
-func iceil(f float64) int {
-	return int(math.Ceil(f))
+	var ix0, iy0, iz0 int
+	ix0, nxSub = ms1.GridSubdomain(domain.Min.X, domain.Max.X, nxDomain, subdomain.Min.X, subdomain.Max.X)
+	iy0, nySub = ms1.GridSubdomain(domain.Min.Y, domain.Max.Y, nyDomain, subdomain.Min.Y, subdomain.Max.Y)
+	iz0, nzSub = ms1.GridSubdomain(domain.Min.Z, domain.Max.Z, nzDomain, subdomain.Min.Z, subdomain.Max.Z)
+	return ix0 + iy0*nxDomain + iz0*(nxDomain+nyDomain), nxSub, nySub, nzSub
 }

md3/grid.go

@@ -158,3 +158,31 @@ func (nra NewtonRaphsonSolver) Root(x0 float32, f func(xGuess float32) float32)
 	}
 	return x_root, -nra.MaxIterations
 }
+
+// GridSubdomain returns the indexing for a 1D subdomain contained within a larger gridded domain of equally spaced nGridPoints points.
+//   - iStart contains the starting index of the first point contained within the subdomain.
+//   - nSub is the number of consecutive grid points contained within the subdomain.
+//
+// This function is primarily used for higher dimension GridSubdomain calls in 2D and 3D packages
+func GridSubdomain(domMin, domMax float32, nGridPoints int, subMin, subMax float32) (iStart, nSub int) {
+	if domMin > domMax || subMin > subMax || subMax > domMax || nGridPoints <= 1 {
+		panic("invalid arguments to gridSubdomain")
+	}
+	dx := (domMax - domMin) / float32(nGridPoints-1) // The size of grid spaces.
+	startOffset := subMin - domMin
+	iStartf := startOffset / dx
+	iStart = min(int(math.Ceil(iStartf)), nGridPoints-1)
+
+	endOffset := subMax - domMin
+	iEndf := endOffset / dx
+	iEnd := min(int(math.Ceil(iEndf)), nGridPoints-1)
+	nSub = iEnd - iStart
+	return iStart, nSub
+}
+
+func min(a, b int) int {
+	if a < b {
+		return a
+	}
+	return b
+}

ms1/ms1.go

@@ -0,0 +1,56 @@
+package ms1
+
+import (
+	"math/rand"
+	"testing"
+)
+
+func Test_gridSubdomainInternal(t *testing.T) {
+	rng := rand.New(rand.NewSource(1))
+	for i := 0; i < 100000; i++ {
+		dmin := float32(rng.Float64())
+		dmax := dmin + float32(rng.Float64())
+		dsize := dmax - dmin
+
+		smin := dmin + float32(rng.Float64())*dsize
+		smaxsize := dmax - smin
+		smax := smin + float32(rng.Float64())*smaxsize
+		if i%8 == 0 {
+			smax = dmax
+		}
+		if i%16 == 0 {
+			smin = dmin
+		}
+		ssize := smax - smin
+		for Nd := 2; Nd < 10; Nd++ {
+			istart, nsub := GridSubdomain(dmin, dmax, Nd, smin, smax)
+			if nsub > Nd {
+				t.Error("subdomain contains more grid points than actual domain", nsub, Nd)
+			} else if istart >= Nd {
+				t.Error("istart >= Nd", istart, Nd)
+			} else if istart+nsub > Nd {
+				t.Error("istart+nsub>Nd", istart, nsub, istart+nsub, Nd)
+			} else if nsub < 0 {
+				t.Error("nsub<0", nsub)
+			}
+			dx := dsize / float32(Nd-1)
+			calc_xmax := dmin + float32(istart+nsub-1)*dx
+			calc_xmin := dmin + float32(istart)*dx
+			if nsub > 0 && calc_xmax > smax {
+				t.Error("xmax > smax", calc_xmax, smax)
+			}
+			if nsub > 0 && calc_xmin < smin {
+				t.Error("xmin < smin", calc_xmin, smin)
+			}
+			if nsub > 0 && calc_xmax > dmax {
+				t.Error("xmax > dmax", calc_xmax, dmax)
+			}
+			_ = ssize
+			if t.Failed() {
+				istart, nsub = GridSubdomain(dmin, dmax, Nd, smin, smax) // For debugging.
+				t.Fatalf("failed on output: i=%d ns=%d  input: d=(%.3f,%.3f) N=%d sub=(%.3f, %.3f) ix=%.3f",
+					istart, nsub, dmin, dmax, Nd, smin, smax, calc_xmin)
+			}
+		}
+	}
+}