Adding ability to specify the inflow velocity profile type

.gitignore

@@ -4,4 +4,8 @@
 *.pyc
 *.egg-info/
 
+*.xdmf
+*.h5
+
 docs/_build/
+results/

Gmsh_code/demo/demo_elipsoidal_shape_concave_surface.txt

@@ -0,0 +1,111 @@
+//+
+SetFactory("OpenCASCADE");
+
+/* create whole box domain */
+Box(1) = {0, 0, 0, 50, 100, 50};
+
+/* create surface on ground (shape shaded by the hill) */
+
+nx = 10;
+ny = 10;
+
+// create all points at front line
+
+p_f = newp;
+For I In {0:nx}
+x = 10+I*30/nx;
+y = 30+(20^2-20^2/15^2*(x-25)^2)^0.5;
+z = 0;
+Point(p_f+I) = {x,y,z,1};
+EndFor
+
+// create all points at rear line
+
+p_b = newp;
+For I In {0:nx-2}
+x = 10+(I+1)*30/nx;
+y = 30-(20^2-20^2/15^2*(x-25)^2)^0.5;
+z = 0;
+Point(p_b+I) = {x,y,z,1};
+EndFor
+
+// create curves at front of bottom surface
+
+l_b_f = newc;
+For I In {0:nx-1}
+Line(l_b_f+I) = {p_f+I, p_f+I+1};
+EndFor
+
+// create curves at back of bottom surface
+
+l_b_b = newc;
+Line(l_b_b) = {p_f, p_b};       
+For I In {1:nx-2}
+Line(l_b_b+I) = {p_b+I-1, p_b+I};
+EndFor
+Line(l_b_b+nx-1) = {p_b+nx-2, p_f+nx};
+
+
+// create curve loop
+
+cl_b = newcl;
+Curve Loop(cl_b) = {l_b_f:l_b_f+nx-1, -l_b_b-nx+1:-l_b_b};
+// create bottom surface
+s_b = news;
+Surface(s_b) = {cl_b};
+
+/* create mountain surface */
+
+p_s = newp;
+For I In {1:nx-1}
+    x = 10+I*30/nx;
+    For J In {1:ny-1}
+
+        // create point for each x
+
+        y = 30-(20^2-20^2/15^2*(x-25)^2)^0.5+2*(20^2-20^2/15^2*(x-25)^2)^0.5/ny*J;
+        z = 20*(1-(x-25)^2/15^2-(y-30)^2/20^2)^0.5;
+        Point(p_s+(ny-1)*(I-1)+J-1) = {x,y,z,1};
+    EndFor
+EndFor
+
+l_s = newc;
+
+For I In {1:nx-1}
+    For J In {1:ny-1}
+        If (J == 1)
+           BSpline(l_s+ny*(I-1)) = {p_b+I-1, p_s+(ny-1)*(I-1)+J-1};
+        Else
+           BSpline(l_s+ny*(I-1)+J-1) = {p_s+(ny-1)*(I-1)+J-2, p_s+(ny-1)*(I-1)+J-1};
+        EndIf
+    EndFor
+    BSpline(l_s+ny*(I-1)+ny-1) = {p_s+(ny-1)*(I-1)+ny-2, p_f+I};
+EndFor
+
+cl_s = newcl;
+For I In {1:nx-1}
+    If (I == 1)
+       Curve Loop(cl_s) = {-l_b_b, l_b_f, -l_s-ny*(I-1)-ny+1:-l_s-ny*(I-1)};
+    Else
+       Curve Loop(cl_s+I-1) = {l_s+ny*(I-2):l_s+ny*(I-2)+ny-1, l_b_f+I-1, -(l_s+ny*(I-1)+ny-1):-(l_s+ny*(I-1)), -l_b_b-I+1};
+    EndIf
+EndFor
+
+Curve Loop(cl_s+nx-1) = {l_b_f+nx-1, -l_b_b-nx+1,  l_s+ny*(nx-2):l_s+ny*(nx-2)+ny-1};
+
+s_s = news;
+For I In {0:nx-1}
+    Printf('%g', I);
+    Surface(s_s+I) = {cl_s+I};
+EndFor
+
+/* create mountain volume */
+
+sl_s = newsl;
+Surface Loop(sl_s) = {s_b, s_s:s_s+nx-1};
+v_m = newv;
+Volume(v_m) = {sl_s};
+
+/* subtract the mountain volume from the whole domain */
+
+BooleanDifference{ Volume{1}; Delete; }{ Volume{v_m}; Delete; }

Gmsh_code/demo/demo_elipsoidal_shape_triangle_surface.geo

@@ -0,0 +1,165 @@
+//+
+SetFactory("OpenCASCADE");
+
+/* Define parameters */
+
+nx = 20;
+ny = 20;
+
+x_min = 10;
+x_max = 40;
+
+y_center = 30;
+
+/* create points on mountain */
+
+p_left = newp;
+Point(p_left) = {x_min, y_center, 0};
+
+p_center = newp;
+
+For I In {1:nx-1}
+    For J In {0:ny}
+        x = x_min+(x_max-x_min)/nx*I;
+        y = 30-(20^2-20^2/15^2*(x-25)^2)^0.5+2*(20^2-20^2/15^2*(x-25)^2)^0.5/ny*J;
+        If (1-(x-25)^2/15^2-(y-30)^2/20^2<0)
+            z = 0;
+        Else
+            z = 20*(1-(x-25)^2/15^2-(y-30)^2/20^2)^0.5;
+        EndIf
+        Point(p_center+(ny+1)*(I-1)+J) = {x, y, z};
+    EndFor
+EndFor
+
+p_right = newp;
+Point(p_right) = {x_max, y_center, 0};
+
+/* Create surface of mountain */
+
+// create lines along y-axis at different x locations
+
+line_ver = newc;
+
+For I In {1:nx-1}
+    For J In {0:ny-1}
+        Line(line_ver+ny*(I-1)+J) = {p_center+(ny+1)*(I-1)+J, p_center+(ny+1)*(I-1)+J+1};
+    EndFor
+EndFor
+
+// create lines at most left surface (connect p_left with p_center+0....p_center+ny)
+
+line_left = newc;
+For J In {0:ny}
+    Line(line_left+J) = {p_left, p_center+J};
+EndFor
+
+// create curve loops at most left surface
+
+cl_left = newcl;
+For J In {0:ny-1}
+    Curve Loop(cl_left+J) = {line_left+J, line_ver+J, -line_left-J-1};
+EndFor
+
+// create most left surfaces
+
+s_left = news;
+For J In {0:ny-1}
+    Surface(s_left+J) = {cl_left+J};
+EndFor
+
+// create lines at center surfaces
+
+line_center = newc;
+For I In {1:nx-2}
+    For J In {0:2*ny}
+        If (J%2==0)
+           Line(line_center+(2*ny+1)*(I-1)+J) = {p_center+(ny+1)*(I-1)+J/2, p_center+(ny+1)*(I)+J/2};
+        Else
+           Line(line_center+(2*ny+1)*(I-1)+J) = {p_center+(ny+1)*(I-1)+(J+1)/2, p_center+(ny+1)*(I)+(J-1)/2};
+        EndIf
+    EndFor
+EndFor
+
+// create line loop at center surfaces
+
+cl_center_left = newcl;
+For I In {1:nx-2}
+    For J In {0:ny-1}
+        Curve Loop(cl_center_left+ny*(I-1)+J) = {line_ver+ny*(I-1)+J, line_center+(2*ny+1)*(I-1)+2*J+1, -(line_center+(2*ny+1)*(I-1)+2*J)};
+    EndFor
+EndFor
+
+s_center_left = news;
+For I In {1:nx-2}
+    For J In {0:ny-1}
+        Surface(s_center_left+ny*(I-1)+J) = {cl_center_left+ny*(I-1)+J};
+    EndFor
+EndFor
+
+
+cl_center_right = newcl;
+For I In {1:nx-2}
+    For J In {0:ny-1}
+        Curve Loop(cl_center_right+ny*(I-1)+J) = {line_center+(2*ny+1)*(I-1)+2*J+1, line_ver+ny*(I)+J, -(line_center+(2*ny+1)*(I-1)+2*J+2)};
+    EndFor
+EndFor
+
+s_center_right = news;
+For I In {1:nx-2}
+    For J In {0:ny-1}
+        Surface(s_center_right+ny*(I-1)+J) = {cl_center_right+ny*(I-1)+J};
+    EndFor
+EndFor
+
+// create lines at most right surface
+
+line_right = newc;
+For J In {0:ny}
+    Line(line_right+J) = {p_right, p_center+(ny+1)*(nx-2)+J};
+EndFor
+
+// create curve loops at most right surface
+
+cl_right = newcl;
+For J In {0:ny-1}
+    Curve Loop(cl_right+J) = {line_right+J, line_ver+ny*(nx-2)+J, -line_right-J-1};
+EndFor
+
+// create most right surfaces
+
+s_right = news;
+For J In {0:ny-1}
+    Surface(s_right+J) = {cl_right+J};
+EndFor
+
+// create bottom surface
+/*
+Printf('%g', line_left);
+Printf('%g', line_center);
+Printf('%g', line_center+(2*ny+1)*(nx-3));
+Printf('%g', -line_right);
+Printf('%g', line_right+ny);
+Printf('%g', -(line_center+(2*ny+1)*(nx-3)+2*ny));
+Printf('%g', -(line_center+2*ny));
+Printf('%g', -line_left-ny);
+*/
+
+
+cl_bottom = newcl;
+Curve Loop(cl_bottom) = {line_left,  line_center:line_center+(2*ny+1)*(nx-3):2*ny+1, -line_right, line_right+ny, -(line_center+(2*ny+1)*(nx-3)+2*ny):-(line_center+2*ny):2*ny+1, -line_left-ny};
+
+s_bottom = news;
+Surface(s_bottom) = {cl_bottom};
+
+// create mountain volume
+
+sc_mountain = news;
+Surface Loop(sc_mountain) = {s_bottom, s_left:s_left+ny-1, s_center_left:s_center_left+ny*(nx-2-1)+ny-1, s_center_right:s_center_right+ny*(nx-2-1)+ny-1, s_right:s_right+ny-1};
+v_m = newv;
+Volume(v_m) = {sc_mountain};
+
+/* create whole box domain */
+v_m_box = newv;
+Box(v_m_box) = {0, 0, 0, 50, 100, 50};
+
+BooleanDifference{ Volume{v_m_box}; Delete; }{ Volume{v_m}; Delete; }