Skip to content

Commit 08e295e

Browse files
IngeborgGjerdeIngeborgGjerde
committed
Cool geometry in coupled 1d3d demo
1 parent c42d7aa commit 08e295e

File tree

1 file changed

+202
-32
lines changed

1 file changed

+202
-32
lines changed

‎demo/Coupled 1d-3d.ipynb‎

Lines changed: 202 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,21 @@
22
"cells": [
33
{
44
"cell_type": "code",
5-
"execution_count": 6,
5+
"execution_count": 1,
66
"metadata": {},
7-
"outputs": [],
7+
"outputs": [
8+
{
9+
"name": "stdout",
10+
"output_type": "stream",
11+
"text": [
12+
"Missing HsMG for fract norm computing\n"
13+
]
14+
}
15+
],
816
"source": [
917
"from fenics import *\n",
1018
"import sys\n",
1119
"sys.path.append('../data/')\n",
12-
"sys.path.append('../graphnics/')\n",
13-
"sys.path.append('../applications/')\n",
14-
"sys.path.append('../../NetworkGen/')\n",
1520
"import matplotlib.pyplot as plt\n",
1621
"import networkx as nx\n",
1722
"import numpy as np\n",
@@ -33,7 +38,7 @@
3338
"- blood flow in vascularized tissue,\n",
3439
"- water flow through the root network in soil,\n",
3540
"- fluid flow through wells drilled in a reservoir.\n",
36-
"- \n",
41+
"\n",
3742
"Each of these application share a common geometry: A network of flow channels $\\Lambda$ (blood vessels, roots or wells) embedded in a surrounding porous media $\\Omega \\subset \\mathbb{R}^3$ (tissue, soil or reservoir).\n",
3843
"\n",
3944
"## Model equations\n",
@@ -42,7 +47,7 @@
4247
"\n",
4348
"$$\n",
4449
"\\begin{align*}\n",
45-
"- \\Delta u &= \\beta(\\hat{u}-\\bar{u})\\delta_\\Gamma \\quad &&\\text{ in } \\Omega, \\\\\n",
50+
"- \\Delta u &= \\beta(\\hat{u}-\\bar{u})\\delta_\\Lambda \\quad &&\\text{ in } \\Omega, \\\\\n",
4651
"- \\partial_{ss} \\hat{u} &= -\\beta(\\hat{u}-\\bar{u}) \\quad &&\\text{ on } \\Lambda\n",
4752
"\\end{align*}\n",
4853
"$$\n",
@@ -71,7 +76,7 @@
7176
},
7277
{
7378
"cell_type": "code",
74-
"execution_count": 7,
79+
"execution_count": 19,
7580
"metadata": {},
7681
"outputs": [],
7782
"source": [
@@ -87,19 +92,12 @@
8792
"mesh3d = UnitCubeMesh(40, 40, 40)\n",
8893
"\n",
8994
"c = mesh3d.coordinates()\n",
95+
"xc, yc, zc = (np.min(node_coords, axis=0)+np.max(node_coords, axis=0))*0.5\n",
96+
"xl, yl, zl = (np.max(node_coords, axis=0)-np.min(node_coords, axis=0))\n",
9097
"\n",
91-
"xmin, ymin, zmin = np.min(node_coords, axis=0)\n",
92-
"xmax, ymax, zmax = np.max(node_coords, axis=0)\n",
93-
"\n",
94-
"c[:, 0] *= (xmax - xmin) * 1.2\n",
95-
"c[:, 1] *= (ymax - ymin) * 1.2\n",
96-
"c[:, 2] *= 0.5\n",
97-
"c[:, 2] -= 0.25\n",
98-
"c[:, 0] -= np.abs(xmin)*1.2\n",
99-
"c[:, 1] -= 0.1\n",
100-
"\n",
101-
"File(\"plots/coupled1d3d/box.pvd\") << Function(FunctionSpace(mesh3d, \"CG\", 1))\n",
102-
"File(\"plots/coupled1d3d/network.pvd\") << Function(FunctionSpace(mesh1d, \"CG\", 1))"
98+
"c[:,:]-=[xc,yc,zc]\n",
99+
"for i, length in enumerate([xl, yl, zl]):\n",
100+
" c[:,i]*=length*1.1"
103101
]
104102
},
105103
{
@@ -113,7 +111,7 @@
113111
},
114112
{
115113
"cell_type": "code",
116-
"execution_count": 8,
114+
"execution_count": 3,
117115
"metadata": {},
118116
"outputs": [],
119117
"source": [
@@ -173,7 +171,7 @@
173171
},
174172
{
175173
"cell_type": "code",
176-
"execution_count": 10,
174+
"execution_count": 4,
177175
"metadata": {},
178176
"outputs": [],
179177
"source": [
@@ -216,15 +214,14 @@
216214
},
217215
{
218216
"cell_type": "code",
219-
"execution_count": 11,
217+
"execution_count": 5,
220218
"metadata": {},
221219
"outputs": [
222220
{
223-
"name": "stdout",
221+
"name": "stderr",
224222
"output_type": "stream",
225223
"text": [
226-
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
227-
"Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
224+
"Averaging over 48 cells: 48it [00:00, 1315.14it/s]\n"
228225
]
229226
}
230227
],
@@ -233,10 +230,8 @@
233230
"a = [[a00, a01], [a10, a11]]\n",
234231
"L = [L0, L1]\n",
235232
"\n",
236-
"W_bcs = [\n",
237-
" [DirichletBC(V3, u_bc_3, \"on_boundary\")],\n",
238-
" [DirichletBC(V1, u_bc_1, \"on_boundary\")],\n",
239-
"]\n",
233+
"W_bcs = [[DirichletBC(V3, u_bc_3, \"on_boundary\")],\n",
234+
" [DirichletBC(V1, u_bc_1, \"on_boundary\")]]\n",
240235
"\n",
241236
"A, b = map(ii_assemble, (a, L))\n",
242237
"A, b = apply_bc(A, b, W_bcs)\n",
@@ -258,13 +253,188 @@
258253
"\n",
259254
"Finally, we can plot the results using paraview. If you open the files `pressure1d.pvd` and `pressure3d.pvd` located in the `plots/coupled1d3d/` folder you should see something like this:\n",
260255
"\n",
261-
"<img alt=\"alt_text\" width=\"500px\" src=\"coupled1d3d.png\" />"
256+
"<img alt=\"alt_text\" width=\"500px\" src=\"coupled1d3d.png\" />\n",
257+
"\n",
258+
"## Coupled 1d-3d flow model for carcinoma\n",
259+
"\n",
260+
"As you might have noticed from the above figure, `fenics_ii` can assemble and solve the coupled 1d-3d flow model even if the 1d mesh does not conform to the 3d mesh. Thus, using e.g. `graphnics` to make the 1d mesh, it is straightforward to extend our approach to more complicated networks. We illustrate this by using a graph constructed from the vascular network of a carcinoma in a rat."
261+
]
262+
},
263+
{
264+
"cell_type": "code",
265+
"execution_count": 48,
266+
"metadata": {},
267+
"outputs": [
268+
{
269+
"name": "stdout",
270+
"output_type": "stream",
271+
"text": [
272+
"Downoading data from https://physiology.arizona.edu/sites/default/files/rattum98_0.txt \n",
273+
"\n",
274+
"Please visit https://physiology.arizona.edu/people/secomb/network for information on how to cite the source for this data\n",
275+
" \n"
276+
]
277+
}
278+
],
279+
"source": [
280+
"# Get 1d network\n",
281+
"from data.secomb_tumours import read_network\n",
282+
"G = read_network.get_graph('https://physiology.arizona.edu/sites/default/files/rattm93b.txt')\n",
283+
"G.make_mesh(4)\n",
284+
"mesh1d = G.global_mesh\n",
285+
"\n",
286+
"# Fit box around it\n",
287+
"pos = nx.get_node_attributes(G, \"pos\")\n",
288+
"node_coords = np.asarray(list(pos.values()))\n",
289+
"\n",
290+
"mesh3d = UnitCubeMesh(10, 10, 10)\n",
291+
"\n",
292+
"c = mesh3d.coordinates()\n",
293+
"xc, yc, zc = (np.min(node_coords, axis=0)+np.max(node_coords, axis=0))*0.5\n",
294+
"xl, yl, zl = (np.max(node_coords, axis=0)-np.min(node_coords, axis=0))\n",
295+
"\n",
296+
"c[:,:]-=[xc,yc,zc]\n",
297+
"for i, length in enumerate([xl, yl, zl]):\n",
298+
" c[:,i]*=length*1.1"
299+
]
300+
},
301+
{
302+
"cell_type": "markdown",
303+
"metadata": {},
304+
"source": [
305+
"Assembling and solving the model is performed as shown earlier, with the exception that the averaging radius is taken from the vasculature radius."
306+
]
307+
},
308+
{
309+
"cell_type": "code",
310+
"execution_count": null,
311+
"metadata": {},
312+
"outputs": [],
313+
"source": [
314+
"W = [FunctionSpace(msh, \"CG\", 1) for msh in [mesh3d, mesh1d]]\n",
315+
"\n",
316+
"u3, u1 = list(map(TrialFunction, W))\n",
317+
"v3, v1 = list(map(TestFunction, W))\n",
318+
"\n",
319+
"\n",
320+
"# Averaging radius is taken from radius data\n",
321+
"class AveragingRadius(UserExpression):\n",
322+
" '''\n",
323+
" Vessel radius as a fenics expression\n",
324+
" '''\n",
325+
" def __init__(self, G, **kwargs):\n",
326+
" self.G = G\n",
327+
" super().__init__(**kwargs)\n",
328+
" def eval(self, value, x):\n",
329+
" p = Point(x[0], x[1], x[2])\n",
330+
" tree = BoundingBoxTree()\n",
331+
" tree.build(mesh1d)\n",
332+
" cell = tree.compute_first_entity_collision(p)\n",
333+
"\n",
334+
" edge_ix = self.G.mf[cell]\n",
335+
" edge = list(G.edges())[edge_ix]\n",
336+
" value[0] = self.G.edges()[edge]['radius']\n",
337+
" \n",
338+
"radius = Radius(G, degree=2)\n",
339+
"\n",
340+
"cylinder = Circle(radius=radius, degree=10)\n",
341+
"\n",
342+
"Pi_u, Pi_v = [Average(func, mesh1d, cylinder) for func in [u3, v3]]"
262343
]
344+
},
345+
{
346+
"cell_type": "markdown",
347+
"metadata": {},
348+
"source": [
349+
"From here we assemble the form and solve:"
350+
]
351+
},
352+
{
353+
"cell_type": "code",
354+
"execution_count": 81,
355+
"metadata": {},
356+
"outputs": [
357+
{
358+
"name": "stdout",
359+
"output_type": "stream",
360+
"text": [
361+
"WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.\n"
362+
]
363+
}
364+
],
365+
"source": [
366+
"beta = Constant(0.000001)\n",
367+
"\n",
368+
"dxGamma = Measure(\"dx\", domain=mesh1d)\n",
369+
"\n",
370+
"# blocks\n",
371+
"a00 = inner(grad(u3), grad(v3)) * dx + circum* beta * inner(Pi_u, Pi_v) * dxGamma\n",
372+
"a01 = -beta * circum * inner(u1, Pi_v) * dxGamma\n",
373+
"a10 = -beta * inner(Pi_u, v1) * dxGamma\n",
374+
"a11 = inner(grad(u1), grad(v1)) * dx + beta * inner(u1, v1) * dxGamma\n",
375+
"\n",
376+
"# right-hand side\n",
377+
"L0 = inner(Constant(0), Pi_v) * dxGamma\n",
378+
"L1 = inner(Constant(0), v1) * dxGamma\n",
379+
"\n",
380+
"\n",
381+
"a = [[a00, a01], [a10, a11]]\n",
382+
"L = [L0, L1]\n",
383+
"\n",
384+
"u_bc_1 = Expression('0.1*x[1]', degree=2)\n",
385+
"\n",
386+
"W_bcs = [[DirichletBC(W[0], u_bc_3, \"on_boundary\")],\n",
387+
" [DirichletBC(W[1], u_bc_1, \"on_boundary\")]]\n",
388+
"\n",
389+
"A, b = map(ii_assemble, (a, L))\n",
390+
"A, b = apply_bc(A, b, W_bcs)\n",
391+
"A, b = map(ii_convert, (A, b))\n",
392+
"\n",
393+
"wh = ii_Function(W)\n",
394+
"solver = LUSolver(A, \"mumps\")\n",
395+
"solver.solve(wh.vector(), b)"
396+
]
397+
},
398+
{
399+
"cell_type": "code",
400+
"execution_count": 82,
401+
"metadata": {},
402+
"outputs": [
403+
{
404+
"name": "stderr",
405+
"output_type": "stream",
406+
"text": [
407+
"Averaging over 1648 cells: 1648it [00:03, 538.39it/s]\n"
408+
]
409+
},
410+
{
411+
"data": {
412+
"text/plain": [
413+
"1"
414+
]
415+
},
416+
"execution_count": 82,
417+
"metadata": {},
418+
"output_type": "execute_result"
419+
}
420+
],
421+
"source": [
422+
"write_vtp(G, [(radius_i, 'radius'), (wh[1], 'p1d')], fname='plots/coupled1d3d/pressure1d.vtp')\n",
423+
"wh[0].rename('p3d', '0.0')\n",
424+
"File('plots/coupled1d3d/pressure3d.pvd') << wh[0]"
425+
]
426+
},
427+
{
428+
"cell_type": "code",
429+
"execution_count": 83,
430+
"metadata": {},
431+
"outputs": [],
432+
"source": []
263433
}
264434
],
265435
"metadata": {
266436
"kernelspec": {
267-
"display_name": "Python 3.6.9 64-bit",
437+
"display_name": "Python 3",
268438
"language": "python",
269439
"name": "python3"
270440
},

0 commit comments

Comments
 (0)