|
51 | 51 | "However, in cases where we have a solution we know that should have no approximation error, we know that the solution should\n",
|
52 | 52 | "be produced to machine precision by the program.\n",
|
53 | 53 | "\n",
|
54 |
| - "The first eight lines of the program are importing the different modules required for solving the problem." |
| 54 | + "We start by importing external modules used alongside `dolfinx`." |
55 | 55 | ]
|
56 | 56 | },
|
57 | 57 | {
|
|
60 | 60 | "metadata": {},
|
61 | 61 | "outputs": [],
|
62 | 62 | "source": [
|
63 |
| - "import dolfinx\n", |
64 | 63 | "import numpy\n",
|
65 | 64 | "import ufl\n",
|
66 | 65 | "from mpi4py import MPI"
|
|
70 | 69 | "cell_type": "markdown",
|
71 | 70 | "metadata": {},
|
72 | 71 | "source": [
|
73 |
| - "A major difference between a traditional FEniCS code and a FEniCSx code, is that one is not advised to use the wildcard import.\n", |
| 72 | + "A major difference between a traditional FEniCS code and a FEniCSx code, is that one is not advised to use the wildcard import. We will see this throughout this first example.\n", |
74 | 73 | "## Generating simple meshes\n",
|
75 |
| - "The next step is to define the discrete domain, _the mesh_\n" |
| 74 | + "The next step is to define the discrete domain, _the mesh_. We do this by importing the built-in `dolfinx.generation.UnitSquareMesh` function, that can build a $[0,1]\\times[0,1]$ mesh, consisting of either triangles or quadrilaterals." |
76 | 75 | ]
|
77 | 76 | },
|
78 | 77 | {
|
|
81 | 80 | "metadata": {},
|
82 | 81 | "outputs": [],
|
83 | 82 | "source": [
|
84 |
| - "mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 8, 8, dolfinx.cpp.mesh.CellType.quadrilateral)" |
| 83 | + "from dolfinx.generation import UnitSquareMesh\n", |
| 84 | + "from dolfinx.mesh import CellType\n", |
| 85 | + "mesh = UnitSquareMesh(MPI.COMM_WORLD, 8, 8, CellType.quadrilateral)" |
85 | 86 | ]
|
86 | 87 | },
|
87 | 88 | {
|
|
111 | 112 | "metadata": {},
|
112 | 113 | "outputs": [],
|
113 | 114 | "source": [
|
114 |
| - " V = dolfinx.FunctionSpace(mesh, (\"CG\", 1))" |
| 115 | + "from dolfinx.fem import FunctionSpace\n", |
| 116 | + "V = FunctionSpace(mesh, (\"CG\", 1))" |
115 | 117 | ]
|
116 | 118 | },
|
117 | 119 | {
|
|
124 | 126 | "and hexahedra). For an overview, see:\n",
|
125 | 127 | "*FIXME: Add link to all the elements we support*\n",
|
126 | 128 | "\n",
|
127 |
| - "The element degree in the code is 1. This means that we are choosing the standard $P_1$ linear Lagrange element, which has degrees of freedom at the vertices. The computed solution will be continuous across elements and linearly varying in $x$ and $y$ inside each element. Higher degree polynomial approximations are obtained by increasing the degree argument. \n", |
| 129 | + "The element degree in the code is 1. This means that we are choosing the standard $P_1$ linear Lagrange element, which has degrees of freedom at the vertices. \n", |
| 130 | + "The computed solution will be continuous across elements and linearly varying in $x$ and $y$ inside each element. Higher degree polynomial approximations are obtained by increasing the degree argument. \n", |
128 | 131 | "\n",
|
129 | 132 | "## Defining the boundary conditions\n",
|
130 | 133 | "\n",
|
131 |
| - "The next step is to specify the boundary condition $u=u_D$ on $\\partial\\Omega_D$, which is done by over several steps. The first step is to define the function $u_D$. Into this function, we would like to interpolate the boundary condition $1 + x^2+2y^2$. We do this by first defining a `dolfinx.Function`, and then using a [lambda-function](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python to define the \n", |
132 |
| - "spatially varying function. As we would like this program to work on one or multiple processors, we have to update the coefficients of $u_D$ that data shared between the processors. We do this by sending (scattering) the ghost values in the underlying data structure of `uD`.\n" |
| 134 | + "The next step is to specify the boundary condition $u=u_D$ on $\\partial\\Omega_D$, which is done by over several steps. \n", |
| 135 | + "The first step is to define the function $u_D$. Into this function, we would like to interpolate the boundary condition $1 + x^2+2y^2$.\n", |
| 136 | + "We do this by first defining a `dolfinx.Function`, and then using a [lambda-function](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) in Python to define the \n", |
| 137 | + "spatially varying function.\n" |
133 | 138 | ]
|
134 | 139 | },
|
135 | 140 | {
|
|
138 | 143 | "metadata": {},
|
139 | 144 | "outputs": [],
|
140 | 145 | "source": [
|
141 |
| - "uD = dolfinx.Function(V)\n", |
142 |
| - "uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)\n", |
143 |
| - "uD.x.scatter_forward()" |
| 146 | + "from dolfinx.fem import Function\n", |
| 147 | + "uD = Function(V)\n", |
| 148 | + "uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)" |
144 | 149 | ]
|
145 | 150 | },
|
146 | 151 | {
|
|
150 | 155 | "We now have the boundary data (and in this case the solution of \n",
|
151 | 156 | "the finite element problem) represented in the discrete function space.\n",
|
152 | 157 | "Next we would like to apply the boundary values to all degrees of freedom that are on the boundary of the discrete domain. We start by identifying the facets (line-segments) representing the outer boundary, using `dolfinx.cpp.mesh.compute_boundary_facets`.\n",
|
153 |
| - "This function returns an array of booleans of the same size as the number of facets on this processor, where `True` indicates that the local facet $i$ is on the boundary. To reduce this to only the indices that are `True`, we use `numpy.where`.\n" |
| 158 | + "This function returns an array of booleans of the same size as the number of facets on this processor, where `True` indicates that the local facet $i$ is on the boundary. To reduce this to only the indices that are `True`, we use `numpy.where`." |
154 | 159 | ]
|
155 | 160 | },
|
156 | 161 | {
|
|
161 | 166 | "source": [
|
162 | 167 | "fdim = mesh.topology.dim - 1\n",
|
163 | 168 | "# Create facet to cell connectivity required to determine boundary facets\n",
|
| 169 | + "from dolfinx.cpp.mesh import compute_boundary_facets\n", |
164 | 170 | "mesh.topology.create_connectivity(fdim, mesh.topology.dim)\n",
|
165 |
| - "boundary_facets = numpy.where(numpy.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]" |
| 171 | + "boundary_facets = numpy.where(numpy.array(compute_boundary_facets(mesh.topology)) == 1)[0]" |
166 | 172 | ]
|
167 | 173 | },
|
168 | 174 | {
|
|
185 | 191 | "metadata": {},
|
186 | 192 | "outputs": [],
|
187 | 193 | "source": [
|
188 |
| - "boundary_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets)\n", |
189 |
| - "bc = dolfinx.DirichletBC(uD, boundary_dofs)" |
| 194 | + "from dolfinx.fem import locate_dofs_topological, DirichletBC\n", |
| 195 | + "boundary_dofs = locate_dofs_topological(V, fdim, boundary_facets)\n", |
| 196 | + "bc = DirichletBC(uD, boundary_dofs)" |
190 | 197 | ]
|
191 | 198 | },
|
192 | 199 | {
|
|
223 | 230 | "metadata": {},
|
224 | 231 | "outputs": [],
|
225 | 232 | "source": [
|
226 |
| - "f = dolfinx.Constant(mesh, -6)" |
| 233 | + "from dolfinx.fem import Constant\n", |
| 234 | + "f = Constant(mesh, -6)" |
227 | 235 | ]
|
228 | 236 | },
|
229 | 237 | {
|
|
273 | 281 | "metadata": {},
|
274 | 282 | "outputs": [],
|
275 | 283 | "source": [
|
276 |
| - "problem = dolfinx.fem.LinearProblem(a, L, bcs=[bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", |
| 284 | + "from dolfinx.fem import LinearProblem\n", |
| 285 | + "problem = LinearProblem(a, L, bcs=[bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", |
277 | 286 | "uh = problem.solve()"
|
278 | 287 | ]
|
279 | 288 | },
|
|
292 | 301 | "metadata": {},
|
293 | 302 | "outputs": [],
|
294 | 303 | "source": [
|
295 |
| - "V2 = dolfinx.FunctionSpace(mesh, (\"CG\", 2))\n", |
296 |
| - "uex = dolfinx.Function(V2)\n", |
297 |
| - "uex.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)\n", |
298 |
| - "uex.x.scatter_forward()" |
| 304 | + "V2 = FunctionSpace(mesh, (\"CG\", 2))\n", |
| 305 | + "uex = Function(V2)\n", |
| 306 | + "uex.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)" |
299 | 307 | ]
|
300 | 308 | },
|
301 | 309 | {
|
302 | 310 | "cell_type": "markdown",
|
303 | 311 | "metadata": {},
|
304 | 312 | "source": [
|
305 |
| - "We compute the error in two different ways. First, we compute the $L^2$-norm of the error, defined by $E=\\sqrt{\\int_\\Omega (u_D-u_h)^2\\mathrm{d} x}$. We use UFL to express the $L^2$-error:" |
| 313 | + "We compute the error in two different ways. First, we compute the $L^2$-norm of the error, defined by $E=\\sqrt{\\int_\\Omega (u_D-u_h)^2\\mathrm{d} x}$. We use UFL to express the $L^2$-error, and use `dolfinx.fem.assemble_scalar` to compute the scalar value." |
306 | 314 | ]
|
307 | 315 | },
|
308 | 316 | {
|
309 | 317 | "cell_type": "code",
|
310 |
| - "execution_count": 12, |
| 318 | + "execution_count": 15, |
311 | 319 | "metadata": {},
|
312 | 320 | "outputs": [],
|
313 | 321 | "source": [
|
| 322 | + "from dolfinx.fem import assemble_scalar\n", |
314 | 323 | "L2_error = ufl.inner(uh - uex, uh - uex) * ufl.dx\n",
|
315 |
| - "error_L2 = numpy.sqrt(dolfinx.fem.assemble_scalar(L2_error))" |
| 324 | + "error_L2 = numpy.sqrt(assemble_scalar(L2_error))" |
316 | 325 | ]
|
317 | 326 | },
|
318 | 327 | {
|
|
329 | 338 | },
|
330 | 339 | {
|
331 | 340 | "cell_type": "code",
|
332 |
| - "execution_count": 13, |
| 341 | + "execution_count": 16, |
333 | 342 | "metadata": {},
|
334 | 343 | "outputs": [
|
335 | 344 | {
|
336 | 345 | "name": "stdout",
|
337 | 346 | "output_type": "stream",
|
338 | 347 | "text": [
|
339 | 348 | "Error_L2 : 8.24e-03\n",
|
340 |
| - "Error_max : 2.22e-15\n" |
| 349 | + "Error_max : 3.11e-15\n" |
341 | 350 | ]
|
342 | 351 | }
|
343 | 352 | ],
|
|
361 | 370 | },
|
362 | 371 | {
|
363 | 372 | "cell_type": "code",
|
364 |
| - "execution_count": 14, |
| 373 | + "execution_count": 17, |
365 | 374 | "metadata": {},
|
366 | 375 | "outputs": [],
|
367 | 376 | "source": [
|
|
371 | 380 | },
|
372 | 381 | {
|
373 | 382 | "cell_type": "code",
|
374 |
| - "execution_count": 15, |
| 383 | + "execution_count": 18, |
375 | 384 | "metadata": {},
|
376 | 385 | "outputs": [],
|
377 | 386 | "source": [
|
|
388 | 397 | },
|
389 | 398 | {
|
390 | 399 | "cell_type": "code",
|
391 |
| - "execution_count": 16, |
| 400 | + "execution_count": 19, |
392 | 401 | "metadata": {},
|
393 | 402 | "outputs": [],
|
394 | 403 | "source": [
|
|
409 | 418 | },
|
410 | 419 | {
|
411 | 420 | "cell_type": "code",
|
412 |
| - "execution_count": 17, |
| 421 | + "execution_count": 20, |
413 | 422 | "metadata": {},
|
414 | 423 | "outputs": [
|
415 | 424 | {
|
416 | 425 | "data": {
|
417 | 426 | "application/vnd.jupyter.widget-view+json": {
|
418 |
| - "model_id": "e67d55e5b10847bc8134ae036945a0bf", |
| 427 | + "model_id": "c8d829e8954c46699aca8c42cd8b5b12", |
419 | 428 | "version_major": 2,
|
420 | 429 | "version_minor": 0
|
421 | 430 | },
|
|
442 | 451 | "cell_type": "markdown",
|
443 | 452 | "metadata": {},
|
444 | 453 | "source": [
|
445 |
| - "## Interactive plotting in notebooks\n", |
446 |
| - "To create an interactive plot using pyvista in a Jupyter notebook we us `pyvista.ITKPlotter` which uses `itkwidgets`. To modify the visualization, click on the three bars in the top left corner." |
| 454 | + "## Alternative plotting\n", |
| 455 | + "To create a more interactive plot using pyvista in a Jupyter notebook we us `pyvista.ITKPlotter` which uses `itkwidgets`. To modify the visualization, click on the three bars in the top left corner." |
447 | 456 | ]
|
448 | 457 | },
|
449 | 458 | {
|
450 | 459 | "cell_type": "code",
|
451 |
| - "execution_count": 18, |
| 460 | + "execution_count": 21, |
452 | 461 | "metadata": {},
|
453 | 462 | "outputs": [
|
454 | 463 | {
|
455 | 464 | "data": {
|
456 | 465 | "application/vnd.jupyter.widget-view+json": {
|
457 |
| - "model_id": "f352aee90a2c403cadaa4c4885907023", |
| 466 | + "model_id": "5a5bdbffdc324a548691f67bb09a5a4c", |
458 | 467 | "version_major": 2,
|
459 | 468 | "version_minor": 0
|
460 | 469 | },
|
|
484 | 493 | },
|
485 | 494 | {
|
486 | 495 | "cell_type": "code",
|
487 |
| - "execution_count": 19, |
| 496 | + "execution_count": 22, |
488 | 497 | "metadata": {},
|
489 |
| - "outputs": [ |
490 |
| - { |
491 |
| - "name": "stderr", |
492 |
| - "output_type": "stream", |
493 |
| - "text": [ |
494 |
| - "2021-09-10 08:49:40.089 ( 4.500s) [main ] VTKFile.cpp:617 WARN| Output data is interpolated into a first order Lagrange space.\n" |
495 |
| - ] |
496 |
| - } |
497 |
| - ], |
| 498 | + "outputs": [], |
498 | 499 | "source": [
|
499 | 500 | "import dolfinx.io\n",
|
500 | 501 | "with dolfinx.io.VTKFile(MPI.COMM_WORLD, \"output.pvd\", \"w\") as vtk:\n",
|
|
537 | 538 | "name": "python",
|
538 | 539 | "nbconvert_exporter": "python",
|
539 | 540 | "pygments_lexer": "ipython3",
|
540 |
| - "version": "3.9.5" |
| 541 | + "version": "3.9.7" |
541 | 542 | }
|
542 | 543 | },
|
543 | 544 | "nbformat": 4,
|
|
0 commit comments