Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 17 additions & 8 deletions chapter2/amr.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -191,6 +191,7 @@
},
"outputs": [],
"source": [
"plotter = pyvista.Plotter()\n",
"curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh))\n",
"curved_grid.cell_data[\"ct\"] = ct.values\n",
"plotter.add_mesh(\n",
Expand Down Expand Up @@ -275,7 +276,8 @@
" else:\n",
" diag_kwargs = {\"diag\": 0.0}\n",
"\n",
" M = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)\n",
" M = dolfinx.fem.petsc.assemble_matrix(\n",
" dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)\n",
" M.assemble()\n",
"\n",
" # Next, we define the SLEPc Eigenvalue Problem Solver (EPS), and set up to use a shift\n",
Expand Down Expand Up @@ -340,15 +342,19 @@
" eta_squared = dolfinx.fem.Function(W)\n",
" f = dolfinx.fem.Constant(mesh, 1.0)\n",
" h = dolfinx.fem.Function(W)\n",
" h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(len(h.x.array), dtype=np.int32))\n",
" h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(\n",
" len(h.x.array), dtype=np.int32))\n",
" n = ufl.FacetNormal(mesh)\n",
"\n",
" G = ( # compute cellwise error estimator\n",
" ufl.inner(h**2 * (f + ufl.div(ufl.grad(uh_r))) ** 2, w) * ufl.dx\n",
" + ufl.inner(h(\"+\") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w(\"+\")) * ufl.dS\n",
" + ufl.inner(h(\"-\") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w(\"-\")) * ufl.dS\n",
" + ufl.inner(h(\"+\") / 2 * ufl.jump(ufl.grad(uh_r), n)\n",
" ** 2, w(\"+\")) * ufl.dS\n",
" + ufl.inner(h(\"-\") / 2 * ufl.jump(ufl.grad(uh_r), n)\n",
" ** 2, w(\"-\")) * ufl.dS\n",
" )\n",
" dolfinx.fem.petsc.assemble_vector(eta_squared.x.petsc_vec, dolfinx.fem.form(G))\n",
" dolfinx.fem.petsc.assemble_vector(\n",
" eta_squared.x.petsc_vec, dolfinx.fem.form(G))\n",
" eta = dolfinx.fem.Function(W)\n",
" eta.x.array[:] = np.sqrt(eta_squared.x.array[:])\n",
"\n",
Expand Down Expand Up @@ -421,12 +427,14 @@
" uh_sign = np.sign(uh_r_min)\n",
" if np.isclose(uh_sign, 0):\n",
" uh_sign = np.sign(uh_r_max)\n",
" assert not np.isclose(uh_sign, 0), \"uh_r has zero values, cannot determine sign.\"\n",
" assert not np.isclose(\n",
" uh_sign, 0), \"uh_r has zero values, cannot determine sign.\"\n",
" uh_r.x.array[:] *= uh_sign\n",
"\n",
" # Update plot with refined mesh\n",
" grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))\n",
" curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space))\n",
" curved_grid = pyvista.UnstructuredGrid(\n",
" *dolfinx.plot.vtk_mesh(uh_r.function_space))\n",
" curved_grid.point_data[\"u\"] = uh_r.x.array\n",
" curved_grid = curved_grid.tessellate()\n",
" curved_actor = plotter.add_mesh(\n",
Expand Down Expand Up @@ -479,7 +487,8 @@
" )\n",
"\n",
" cells_to_mark = mark_cells(uh_r, lam)\n",
" mesh, (_, ft) = geoModel.refineMarkedElements(mesh.topology.dim, cells_to_mark)\n",
" mesh, (_, ft) = geoModel.refineMarkedElements(\n",
" mesh.topology.dim, cells_to_mark)\n",
" curved_mesh = geoModel.curveField(order)\n",
" write_frame(plotter, uh_r)\n",
"\n",
Expand Down
25 changes: 17 additions & 8 deletions chapter2/amr.py
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,7 @@
# Again, we visualize the curved mesh with pyvista.

# + tags=["hide-input"]
plotter = pyvista.Plotter()
curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(curved_mesh))
curved_grid.cell_data["ct"] = ct.values
plotter.add_mesh(
Expand Down Expand Up @@ -157,7 +158,8 @@ def solve(
else:
diag_kwargs = {"diag": 0.0}

M = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)
M = dolfinx.fem.petsc.assemble_matrix(
dolfinx.fem.form(m), bcs=[bc], **diag_kwargs)
M.assemble()

# Next, we define the SLEPc Eigenvalue Problem Solver (EPS), and set up to use a shift
Expand Down Expand Up @@ -208,15 +210,19 @@ def mark_cells(uh_r: dolfinx.fem.Function, lam: float):
eta_squared = dolfinx.fem.Function(W)
f = dolfinx.fem.Constant(mesh, 1.0)
h = dolfinx.fem.Function(W)
h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(len(h.x.array), dtype=np.int32))
h.x.array[:] = mesh.h(mesh.topology.dim, np.arange(
len(h.x.array), dtype=np.int32))
n = ufl.FacetNormal(mesh)

G = ( # compute cellwise error estimator
ufl.inner(h**2 * (f + ufl.div(ufl.grad(uh_r))) ** 2, w) * ufl.dx
+ ufl.inner(h("+") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("+")) * ufl.dS
+ ufl.inner(h("-") / 2 * ufl.jump(ufl.grad(uh_r), n) ** 2, w("-")) * ufl.dS
+ ufl.inner(h("+") / 2 * ufl.jump(ufl.grad(uh_r), n)
** 2, w("+")) * ufl.dS
+ ufl.inner(h("-") / 2 * ufl.jump(ufl.grad(uh_r), n)
** 2, w("-")) * ufl.dS
)
dolfinx.fem.petsc.assemble_vector(eta_squared.x.petsc_vec, dolfinx.fem.form(G))
dolfinx.fem.petsc.assemble_vector(
eta_squared.x.petsc_vec, dolfinx.fem.form(G))
eta = dolfinx.fem.Function(W)
eta.x.array[:] = np.sqrt(eta_squared.x.array[:])

Expand Down Expand Up @@ -252,12 +258,14 @@ def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function):
uh_sign = np.sign(uh_r_min)
if np.isclose(uh_sign, 0):
uh_sign = np.sign(uh_r_max)
assert not np.isclose(uh_sign, 0), "uh_r has zero values, cannot determine sign."
assert not np.isclose(
uh_sign, 0), "uh_r has zero values, cannot determine sign."
uh_r.x.array[:] *= uh_sign

# Update plot with refined mesh
grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
curved_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(uh_r.function_space))
curved_grid = pyvista.UnstructuredGrid(
*dolfinx.plot.vtk_mesh(uh_r.function_space))
curved_grid.point_data["u"] = uh_r.x.array
curved_grid = curved_grid.tessellate()
curved_actor = plotter.add_mesh(
Expand Down Expand Up @@ -296,7 +304,8 @@ def write_frame(plotter: pyvista.Plotter, uh_r: dolfinx.fem.Function):
)

cells_to_mark = mark_cells(uh_r, lam)
mesh, (_, ft) = geoModel.refineMarkedElements(mesh.topology.dim, cells_to_mark)
mesh, (_, ft) = geoModel.refineMarkedElements(
mesh.topology.dim, cells_to_mark)
curved_mesh = geoModel.curveField(order)
write_frame(plotter, uh_r)

Expand Down
6 changes: 3 additions & 3 deletions docker/Dockerfile
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ RUN curl -sL https://deb.nodesource.com/setup_18.x -o nodesource_setup.sh && \

# Install netgen from source
RUN apt-get update && \
apt-get -y install python3 python3-tk libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev liblapacke-dev libocct-data-exchange-dev libocct-draw-dev occt-misc libtbb-dev libxi-dev
apt-get -y install python3 python3-tk libpython3-dev libxmu-dev tk-dev tcl-dev cmake git g++ libglu1-mesa-dev liblapacke-dev libocct-data-exchange-dev libocct-draw-dev occt-misc libtbb-dev libxi-dev
WORKDIR /ngsuite
RUN git clone --branch=v6.2.2505 --single-branch https://github.com/NGSolve/netgen.git netgen-src
WORKDIR /ngsuite/netgen-src
Expand All @@ -39,9 +39,9 @@ RUN python3 -m pip install vtk
ADD pyproject.toml /tmp/pyproject.toml
WORKDIR /tmp
# As we install netgen from source we don't need the netgen deps here
RUN python3 -m pip install --no-cache-dir --no-binary=h5py -v .
RUN python3 -m pip install --no-cache-dir --no-binary=h5py --no-build-isolation -v .
RUN python3 -m pip cache purge

RUN python3 -m pip install ngspetsc --no-deps --no-build-isolation
ENV PYVISTA_OFF_SCREEN="false"
ENV PYVISTA_JUPYTER_BACKEND="trame"
ENV LIBGL_ALWAYS_SOFTWARE=1
Expand Down