Run Rainier on aniso mesh it good

Author: danshapero

notebooks/rainier-monolithic.ipynb

@@ -12,7 +12,7 @@
     "from numpy import pi as π\n",
     "import matplotlib.pyplot as plt\n",
     "import firedrake\n",
-    "from firedrake import inner, grad, dx, ds, exp, min_value, max_value, Constant, derivative\n",
+    "from firedrake import inner, grad, dx, ds, exp, min_value, max_value, Constant\n",
     "import irksome\n",
     "import icepack\n",
     "from icepack2 import model"
@@ -34,6 +34,48 @@
     "We'll use a computational domain about half that size."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "de2a5adf-a9b1-4948-aa16-5f7459f670c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def read_mesh_file(mesh_file):\n",
+    "    lines = mesh_file.readlines()\n",
+    "    line_number = [\"Vertices\" in line for line in lines].index(True) + 1\n",
+    "    num_vertices = int(lines[line_number])\n",
+    "    points = np.zeros((num_vertices, 2))\n",
+    "    vertex_lines = lines[line_number + 1: line_number + num_vertices + 1]\n",
+    "    for index, line in enumerate(vertex_lines):\n",
+    "        points[index, :] = [float(val) for val in line.split(\" \")[:2]]\n",
+    "\n",
+    "    line_number = [\"Triangles\" in line for line in lines].index(True) + 1\n",
+    "    num_triangles = int(lines[line_number])\n",
+    "    triangles = np.zeros((num_triangles, 3), dtype=int)\n",
+    "    triangle_lines = lines[line_number + 1: line_number + num_triangles + 1]\n",
+    "    for index, line in enumerate(triangle_lines):\n",
+    "        triangles[index, :] = [int(val) for val in line.split(\" \")[:3]]\n",
+    "\n",
+    "    return points, triangles - 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f071b37-552a-4096-b36f-5782d6d53dc9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mesh_filename = \"disk_anisotropic.mesh\"\n",
+    "with open(mesh_filename, \"r\") as mesh_file:\n",
+    "    points, triangles = read_mesh_file(mesh_file)\n",
+    "\n",
+    "comm = firedrake.COMM_WORLD\n",
+    "plex = firedrake.mesh.plex_from_cell_list(2, triangles, points, comm)\n",
+    "mesh = firedrake.Mesh(plex)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -42,10 +84,21 @@
    "outputs": [],
    "source": [
     "radius = Constant(12e3)\n",
-    "mesh = firedrake.UnitDiskMesh(4)\n",
     "mesh.coordinates.dat.data[:] *= float(radius)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1019369b-c328-40be-b6e1-8f275aea433e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "ax.set_aspect(\"equal\")\n",
+    "firedrake.triplot(mesh, axes=ax);"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "55e5e816-7f40-435e-aa22-78051037a76b",