encapsulate boundaries

3 files changed, 109 insertions, 34 deletions

diff --git a/main.c b/main.c
index 75f6446..275fea3 100644
--- a/main.c
+++ b/main.c
@@ -32,6 +32,17 @@ int main (void)
 	double yD = -1.0; // lower boundary
 	double yU =  1.0; // upper boundary
 
+	// boundary type
+	int north = 0;
+	int south = 1;
+	int east  = 2;
+	int west  = 2;
+
+	double val_north = 40.0;
+	double val_south = 1.8e-5;
+	double val_east  = 0.0;
+	double val_west  = 0.0;
+
 	double dx = (xR - xL) / nx;
 	double dy = (yU - yD) / ny;
 
@@ -46,7 +57,6 @@ int main (void)
 	double* a  = (double*) calloc(5 * dim, sizeof(double));
 	double* c  = (double*) calloc(dim, sizeof(double));
 	double* r  = (double*) calloc(dim, sizeof(double));
-	double* ar = (double*) calloc(dim, sizeof(double));
 	double* x  = (double*) calloc(dim, sizeof(double));
 	double* vx  = (double*) calloc(dim, sizeof(double));
 	double* vy  = (double*) calloc(dim, sizeof(double));
@@ -78,39 +88,8 @@ int main (void)
 
 	// initialize boundary conditions
 	//
-
-
-	// dirichlet boundary
-	for (int i = 0; i < nx; i ++) {
-		x[i] = 40.0;
-	}
-
-	multiply(&m, x, ar);
-
-	for (int i = 0; i < dim; i ++) {
-		r[i] = r[i] - ar[i];
-		if (i < nx) {
-			r[i] = x[i];
-		}
-	}
-
-	free(ar);
-
-	for (int i = 0; i < nx; i ++) {
-		for (int j = 0; j < dim; j ++) {
-			if (i != j) {
-				put(&m, j, i, 0.0);
-				put(&m, i, j, 0.0);
-			} else {
-				put(&m, j, i, 1.0);
-			}
-		}
-	}
-
-	// neumann boundary
-	for (int i = dim - nx; i < dim; i ++) {
-		r[i] = 1.8e-5;
-	}
+	boundary(dim, nx, ny, north, south, east, west, val_north, val_south, 
+		val_east, val_west, &m, r, x);
 
 	// solve the system A*x = b
 	//
diff --git a/sysmat.c b/sysmat.c
index d0c172b..d390c37 100644
--- a/sysmat.c
+++ b/sysmat.c
@@ -161,6 +161,99 @@ void system_matrix(struct diag* m, int nx, int ny, double* val)
 
 }
 
+
+/**
+ * set boundary conditions
+ *
+ * inout: m
+ **/
+void boundary(int dim, int nx, int ny, int north, int south, int east, 
+	int west, double val_north, double val_south, double val_east, 
+	double val_west, struct diag* m, double* r, double* x)
+{
+
+	// 0: dirichlet | 1: neumann | 2: closed wall
+
+
+
+	if (north == 0) {
+
+		for (int i = 0; i < nx;  i ++) {
+			x[i] = val_north;
+		}
+
+		double* ar = (double*) calloc(dim, sizeof(double));
+
+		multiply(m, x, ar);
+
+		for (int i = 0; i < dim; i ++) {
+			r[i] = r[i] - ar[i];
+			if (i < nx) {
+				r[i] = x[i];
+			}
+		}
+
+		free(ar);
+
+		for (int i = 0; i < nx; i ++) {
+			for (int j = 0; j < dim; j ++) {
+				if (i != j) {
+					put(m, j, i, 0.0);
+					put(m, i, j, 0.0);
+				} else {
+					put(m, j, i, 1.0);
+				}
+			}
+		}
+
+	} else if (north == 1) {
+
+		for (int i = 0; i < nx; i ++) {
+			r[i] = val_north;
+		}
+
+	}
+
+	if (south == 0) {
+
+		for (int i = dim - nx; i < dim; i ++) {
+			x[i] = val_south;
+		}
+
+		double* ar = (double*) calloc(dim, sizeof(double));
+
+		multiply(m, x, ar);
+
+		for (int i = 0; i < dim; i ++) {
+			r[i] = r[i] - ar[i];
+			if ((i >= dim - nx) && (i < dim)) {
+				r[i] = x[i];
+			}
+		}
+
+		free(ar);
+
+		for (int i = dim - nx; i < dim; i ++) {
+			for (int j = 0; j < dim; j ++) {
+				if (i != j) {
+					put(m, j, i, 0.0);
+					put(m, i, j, 0.0);
+				} else {
+					put(m, j, i, 1.0);
+				}
+			}
+		}
+
+	} else if (south == 1) {
+
+		for (int i = dim - nx; i < dim; i ++) {
+			r[i] = val_south;
+		}
+
+	}
+
+}
+
 /**
  * compute darcy velocities
  *
diff --git a/sysmat.h b/sysmat.h
index 3837d38..bcc2514 100644
--- a/sysmat.h
+++ b/sysmat.h
@@ -16,6 +16,9 @@
 
 void system_matrix(struct diag*, int, int, double*);
 
+void boundary(int, int, int, int, int, int, int, double, double, double, 
+	double, struct diag*, double*, double*);
+
 void darcy(double*, double*, int, int, double, double, double*, double*);
 
 double harmonic(double, double);