Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
29cc8e1
new exercise files
campospinto Mar 11, 2025
186494b
pbm formulation
campospinto Mar 11, 2025
7fcf870
add sheet1 exercise 2.2
jowezarek Mar 24, 2025
02e7cc5
fix import typo
jowezarek Mar 24, 2025
f347534
fix import type 2
jowezarek Mar 24, 2025
bdfac39
minor beauty updates
jowezarek Mar 25, 2025
31a2c12
add sheet1 exercise 1.4
jowezarek Mar 31, 2025
cada87e
modify rhs of sheet1 exercise 1.3 (prev. exercise 2.2)
jowezarek Mar 31, 2025
6dcdb41
add sheet2 exercise 2.3
jowezarek Mar 31, 2025
cb5f634
add sheet2 exercise 2.2
jowezarek Mar 31, 2025
4ab8711
Merge branch 'master' into exercises
jowezarek Mar 31, 2025
e92e5ec
add sheet3 exercise 3.2
jowezarek Apr 1, 2025
4e18765
Merge branch 'exercises' of https://github.com/pyccel/IGA-Python into…
jowezarek Apr 1, 2025
f2f6b20
add (first version of) sheet4 exercise 4.7
jowezarek Apr 2, 2025
f9537c0
complete sheet4 exercise 4.7
jowezarek Apr 3, 2025
763405b
add sheet1 exercise 1.3 iii)
jowezarek Apr 3, 2025
b9b5399
change .vector_space to .coeff_space
jowezarek Apr 4, 2025
15c5a39
minor tidy up
jowezarek Apr 4, 2025
b9f23ce
Merge branch 'master' into exercises
yguclu Jul 16, 2025
6a38e5f
update the exercises corresponding to exercise sheet 1 (2026)
jowezarek Apr 21, 2026
a32b9a5
update the exercises corresponding to exercise sheet 2 (2026)
jowezarek May 5, 2026
b0c603f
update the exercises corresponding to exercise sheet 3
jowezarek May 10, 2026
81aef1f
update the exercises corresponding to exercise sheet 4
jowezarek Jun 5, 2026
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
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -19,3 +19,5 @@ usr
.env
.iga-python
*_autosummary*

.DS_Store
4 changes: 4 additions & 0 deletions _toc.yml
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,10 @@ parts:
- file: chapter3/navier-stokes-steady
sections:
- file: chapter3/navier-stokes-steady-streamfunction-velocity
- caption: Exercises
chapters:
- file: exercises/harmonic-fem
- file: exercises/harmonic-feec
- caption: Advanced subjects
chapters:
- file: chapter4/subdomains
Expand Down
285 changes: 285 additions & 0 deletions exercises/electro_static_2d.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,285 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Note: This notebook must be changed by the student!**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2D Solver for Electro-Static Problem on the Unit Square\n",
"\n",
"In this exercise we write a solver for the 2D electro-static problem on the unit square, which can be shown to be equivalent to an $\\mathcal{L}^1$ Hodge-Laplace problem.\n",
"\n",
"\\begin{align*}\n",
" &\\text{Find }\\boldsymbol{E} \\in D(\\mathcal{L}^1) \\\\\n",
" &\\mathcal{L}^1\\boldsymbol{E} = (\\boldsymbol{curl}\\ curl - \\boldsymbol{grad}\\ div)\\boldsymbol{E} = -\\boldsymbol{grad}\\ \\rho \\quad \\text{ in }\\Omega=(0,1)^2,\n",
"\\end{align*}\n",
"for the specific choice\n",
"\\begin{align*}\n",
" \\rho(x, y) = 2\\sin(\\pi x)\\sin(\\pi y)\n",
"\\end{align*}\n",
"\n",
"## The Discrete Variational Formulation\n",
"\n",
"The corresponding discrete variational formulation reads\n",
"\n",
"\\begin{align*}\n",
" &\\text{Find }\\boldsymbol{E} \\in V_h^1 \\subset V^1 = H_0(curl;\\Omega)\\text{ such that } \\\\\n",
" &a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V_h^1\n",
"\\end{align*}\n",
"\n",
"where \n",
"- $a(\\boldsymbol{E},\\boldsymbol{v}) \\coloneqq \\int_{\\Omega} (curl\\ \\boldsymbol{E})(curl\\ \\boldsymbol{v}) + (\\widetilde{div}_h\\boldsymbol{E})(\\widetilde{div}_h\\boldsymbol{v}) ~ d\\Omega$ ,\n",
"- $l(\\boldsymbol{v}) := -\\int_{\\Omega} \\boldsymbol{grad}\\ \\rho\\cdot\\boldsymbol{v} ~ d\\Omega$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discrete Model using de Rham objects"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from sympde.calculus import inner\n",
"from sympde.topology import elements_of, Square, Derham\n",
"from sympde.expr.expr import BilinearForm, integral\n",
"\n",
"domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n",
"derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n",
"\n",
"V0 = derham.V0\n",
"V1 = derham.V1\n",
"V2 = derham.V2\n",
"\n",
"u0, v0 = elements_of(V0, names='u0, v0')\n",
"u1, v1 = elements_of(V1, names='u1, v1')\n",
"u2, v2 = elements_of(V2, names='u2, v2')\n",
"\n",
"# bilinear forms corresponding to V0, V1 and V2 mass matrices\n",
"m0 = BilinearForm((u0, v0), integral(domain, u0 * v0))\n",
"m1 = BilinearForm((u1, v1), integral(domain, inner(u1, v1)))\n",
"m2 = BilinearForm((u2, v2), integral(domain, u2 * v2))\n",
"\n",
"# callable source term \\rho\n",
"from sympy import pi, sin, cos, lambdify, Tuple, Matrix\n",
"x,y = domain.coordinates\n",
"rho = 2*sin(pi*x)*sin(pi*y)\n",
"rho_lambdified = lambdify((x, y), rho)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from psydac.api.discretization import discretize\n",
"from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n",
"\n",
"backend = PSYDAC_BACKEND_GPYCCEL"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ncells = [32, 32] # Bspline cells\n",
"degree = [4, 4] # Bspline degree"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# discretize domain and derham\n",
"# <-\n",
"# <-\n",
"\n",
"# define FEM spaces V0_h, V1_h and V2_h\n",
"#V0_h = derham_h.V0\n",
"#V1_h = derham_h.V1\n",
"#V2_h = derham_h.V2\n",
"\n",
"# Commuting projection operators\n",
"#P0, P1, P2 = derham_h.projectors()\n",
"\n",
"# Exterior derivative operators (grad and curl)\n",
"# <-\n",
"\n",
"# Mass matrices (discretization and assembly)\n",
"# <-\n",
"# <-\n",
"# <-\n",
"\n",
"# <-\n",
"# <-\n",
"# <-\n",
"\n",
"# Boundary Projectors\n",
"from psydac.fem.projectors import DirichletProjector\n",
"from psydac.linalg.basic import IdentityOperator\n",
"\n",
"#P_H1 = DirichletProjector(V0_h)\n",
"#P_Hcurl = DirichletProjector(V1_h)\n",
"\n",
"#I0 = IdentityOperator(V0_h.coeff_space)\n",
"#I1 = IdentityOperator(V1_h.coeff_space)\n",
"\n",
"#P_H1_Gamma = I0 - P_H1\n",
"#P_Hcurl_Gamma = I1 - P_Hcurl\n",
"\n",
"# The inner product < div v, div E > poses a difficulty for the projection method.\n",
"# We need to construct the \"regularized\" mass matrix M0_0 \n",
"#M0_0 = P_H1 @ M0 @ P_H1 + P_H1_Gamma\n",
"\n",
"# Verify, that in matrix form the weak divergence reads\n",
"# weak_div = - inverse(M0_0) @ P_H1 @ @ Grad.T @ M1 \n",
"# Conclude, that the matrix representation of the bilinear form \n",
"# (u, v) \\mapsto L2-inner product of weak_div(u) and weak_div(v)\n",
"# takes the form\n",
"# M1 @ Grad @ inverse(M0_0) @ P_H1 @ Grad.T @ M1\n",
"\n",
"# Assemble the system matrix\n",
"from psydac.linalg.solvers import inverse\n",
"#M0_0_inv = inverse(M0_0, 'CG', maxiter=1000, tol=1e-11)\n",
"\n",
"# System Matrix A and A_bc\n",
"#A = # <-\n",
"#A_bc = P_Hcurl @ A @ P_Hcurl + P_Hcurl_Gamma\n",
"\n",
"# rhs vector f and f_bc\n",
"#rho_coeffs = P0(rho_lambdified).coeffs\n",
"#f = # <-\n",
"#f_bc = P_Hcurl @ f"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import time\n",
"\n",
"tol = 1e-10\n",
"maxiter = 1000\n",
"\n",
"#A_bc_inv = inverse(A_bc, 'CG', tol=tol, maxiter=maxiter)\n",
"\n",
"t0 = time.time()\n",
"#E_h = A_bc_inv @ f_bc\n",
"t1 = time.time()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing the error norm\n",
"\n",
"As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n",
"In this example, the analytical solution is given by\n",
"\n",
"$$\n",
"\\boldsymbol{E}_{ex}(x, y) = -\\frac{1}{\\pi}\\left(\\begin{matrix} \n",
" \\cos(\\pi x) \\sin(\\pi y) \\\\\n",
" \\cos(\\pi y) \\sin(\\pi x)\n",
" \\end{matrix}\\right)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from psydac.fem.basic import FemField\n",
"\n",
"from sympde.expr import Norm\n",
"\n",
"E_ex = Tuple( (-1/pi)*cos(pi*x) * sin(pi*y),\n",
" (-1/pi)*cos(pi*y) * sin(pi*x) )\n",
"#E_h_FemField = FemField(V1_h, E_h)\n",
"\n",
"error = Matrix([u1[0] - E_ex[0], u1[1] - E_ex[1]])\n",
"\n",
"# create the formal Norm object\n",
"l2norm = Norm(error, domain, kind='l2')\n",
"\n",
"# discretize the norm\n",
"#l2norm_h = discretize(l2norm, domain_h, V1_h, backend=backend)\n",
"\n",
"# assemble the norm\n",
"#l2_error = l2norm_h.assemble(u1=E_h_FemField)\n",
"\n",
"# print the result\n",
"#print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
"#print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n",
"#print( '> CG info :: ',A_bc_inv.get_info() )\n",
"#print( '> L2 error :: {:.2e}'.format( l2_error ) )\n",
"#print( '' )\n",
"#print( '> Solution time :: {:.3g}'.format( t1-t0 ) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from utils import plot\n",
"\n",
"E_ex_x = lambdify((x, y), E_ex[0])\n",
"E_ex_y = lambdify((x, y), E_ex[1])\n",
"#E_h_x = E_h_FemField[0]\n",
"#E_h_y = E_h_FemField[1]\n",
"#error_x = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n",
"#error_y = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n",
"\n",
"#plot(gridsize_x = 100, \n",
"# gridsize_y = 100, \n",
"# title = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n",
"# funs = [E_ex_x, E_h_x, error_x], \n",
"# titles = [r'$(\\boldsymbol{E}_{ex})_1(x,y)$', r'$(\\boldsymbol{E}_h)_1(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_1(x,y)|$'],\n",
"# surface_plot = True\n",
"#)\n",
"\n",
"#plot(gridsize_x = 100,\n",
"# gridsize_y = 100,\n",
"# title = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n",
"# funs = [E_ex_y, E_h_y, error_y],\n",
"# titles = [r'$(\\boldsymbol{E}_{ex})_2(x,y)$', r'$(\\boldsymbol{E}_h)_2(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_2(x,y)|$'],\n",
"# surface_plot = True\n",
"#)"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 4
}
Loading
Loading