Skip to content

Commit 29cc8e1

Browse files
campospintocampospinto
committed
new exercise files
1 parent babf0d4 commit 29cc8e1

File tree

3 files changed

+552
-0
lines changed

3 files changed

+552
-0
lines changed

‎_toc.yml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,10 @@ parts:
6666
- file: chapter3/navier-stokes-steady
6767
sections:
6868
- file: chapter3/navier-stokes-steady-streamfunction-velocity
69+
- caption: Exercises
70+
chapters:
71+
- file: exercises/harmonic-fem
72+
- file: exercises/harmonic-feec
6973
- caption: Advanced subjects
7074
chapters:
7175
- file: chapter4/subdomains

‎exercises/harmonic-feec.ipynb‎

Lines changed: 274 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,274 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "9f28f9af",
6+
"metadata": {},
7+
"source": [
8+
"# Harmonic fields with structure-preserving (feec) fem spaces\n",
9+
"\n",
10+
"In this example we .... consider the vector Poisson equation with homogeneous Dirichlet boundary conditions:\n",
11+
"\n",
12+
"$$\n",
13+
"\\begin{align}\n",
14+
" - \\nabla^2 \\mathbf{u} = \\mathbf{f} \\quad \\mbox{in} ~ \\Omega, \\quad \\quad \n",
15+
" \\mathbf{u} = 0 \\quad \\mbox{on} ~ \\partial \\Omega.\n",
16+
"\\end{align}\n",
17+
"$$\n",
18+
"\n",
19+
"## The Variational Formulation\n",
20+
"\n",
21+
"The corresponding variational formulation, using $\\mathbf{H}^1$ formulation, *i.e.* all components are in $H^1$, reads \n",
22+
"\n",
23+
"$$\n",
24+
"\\begin{align}\n",
25+
" \\text{find $\\mathbf{u} \\in V$ such that} \\quad \n",
26+
" a(\\mathbf{u},\\mathbf{v}) = l(\\mathbf{v}) \\quad \\forall \\mathbf{v} \\in V,\n",
27+
"\\end{align}\n",
28+
"$$\n",
29+
"\n",
30+
"where \n",
31+
"\n",
32+
"- $V \\subset \\mathbf{H}_0^1(\\Omega)$, \n",
33+
"- $a(\\mathbf{u},\\mathbf{v}) := \\int_{\\Omega} \\nabla \\mathbf{u} : \\nabla \\mathbf{v} ~ d\\Omega$,\n",
34+
"- $l(\\mathbf{v}) := \\int_{\\Omega} \\mathbf{f} \\cdot \\mathbf{v} ~ d\\Omega$."
35+
]
36+
},
37+
{
38+
"cell_type": "markdown",
39+
"id": "5a958607",
40+
"metadata": {},
41+
"source": [
42+
"## Formal Model"
43+
]
44+
},
45+
{
46+
"cell_type": "code",
47+
"execution_count": null,
48+
"id": "d742586c",
49+
"metadata": {},
50+
"outputs": [],
51+
"source": [
52+
"from sympde.expr import BilinearForm, LinearForm, integral\n",
53+
"from sympde.expr import find, EssentialBC, Norm, SemiNorm\n",
54+
"from sympde.topology import VectorFunctionSpace, Cube, element_of\n",
55+
"from sympde.calculus import grad, inner, dot\n",
56+
"\n",
57+
"from sympy import pi, sin, Tuple, Matrix\n",
58+
"\n",
59+
"domain = Cube()\n",
60+
"\n",
61+
"V = VectorFunctionSpace('V', domain)\n",
62+
"\n",
63+
"x,y,z = domain.coordinates\n",
64+
"\n",
65+
"u,v = [element_of(V, name=i) for i in ['u', 'v']]\n",
66+
"\n",
67+
"# bilinear form\n",
68+
"a = BilinearForm((u,v), integral(domain , inner(grad(v), grad(u))))\n",
69+
"\n",
70+
"# linear form\n",
71+
"f1 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
72+
"f2 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
73+
"f3 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
74+
"f = Tuple(f1, f2, f3)\n",
75+
"\n",
76+
"l = LinearForm(v, integral(domain, dot(f,v)))\n",
77+
"\n",
78+
"# Dirichlet boundary conditions\n",
79+
"bc = [EssentialBC(u, 0, domain.boundary)]\n",
80+
"\n",
81+
"# Variational problem\n",
82+
"equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)"
83+
]
84+
},
85+
{
86+
"cell_type": "markdown",
87+
"id": "4f983ece",
88+
"metadata": {},
89+
"source": [
90+
"## Discretization"
91+
]
92+
},
93+
{
94+
"cell_type": "markdown",
95+
"id": "51095918",
96+
"metadata": {},
97+
"source": [
98+
"We shall need the **discretize** function from **PsyDAC**."
99+
]
100+
},
101+
{
102+
"cell_type": "code",
103+
"execution_count": null,
104+
"id": "a2a0a2a1",
105+
"metadata": {},
106+
"outputs": [],
107+
"source": [
108+
"from psydac.api.discretization import discretize"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"id": "00e54163",
115+
"metadata": {},
116+
"outputs": [],
117+
"source": [
118+
"degree = [2,2,2]\n",
119+
"ncells = [8,8,8]"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"id": "5999c62b",
126+
"metadata": {},
127+
"outputs": [],
128+
"source": [
129+
"# Create computational domain from topological domain\n",
130+
"domain_h = discretize(domain, ncells=ncells, comm=None)\n",
131+
"\n",
132+
"# Create discrete spline space\n",
133+
"Vh = discretize(V, domain_h, degree=degree)\n",
134+
"\n",
135+
"# Discretize equation\n",
136+
"equation_h = discretize(equation, domain_h, [Vh, Vh])"
137+
]
138+
},
139+
{
140+
"cell_type": "markdown",
141+
"id": "7b29fbcf",
142+
"metadata": {},
143+
"source": [
144+
"## Solving the PDE"
145+
]
146+
},
147+
{
148+
"cell_type": "code",
149+
"execution_count": null,
150+
"id": "541192ee",
151+
"metadata": {},
152+
"outputs": [],
153+
"source": [
154+
"uh = equation_h.solve()"
155+
]
156+
},
157+
{
158+
"cell_type": "markdown",
159+
"id": "5174c4b5",
160+
"metadata": {},
161+
"source": [
162+
"## Computing the error norm\n",
163+
"\n",
164+
"When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n",
165+
"SymPDE allows you to do so, by creating the **Norm** object.\n",
166+
"In this example, the analytical solution is given by\n",
167+
"\n",
168+
"$$\n",
169+
"u_e = \\sin(\\pi x) \\sin(\\pi y) \\sin(\\pi z)\n",
170+
"$$"
171+
]
172+
},
173+
{
174+
"cell_type": "markdown",
175+
"id": "3a31c46f",
176+
"metadata": {},
177+
"source": [
178+
"### Computing the $L^2$ norm"
179+
]
180+
},
181+
{
182+
"cell_type": "code",
183+
"execution_count": null,
184+
"id": "5925c6cd",
185+
"metadata": {},
186+
"outputs": [],
187+
"source": [
188+
"ue1 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
189+
"ue2 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
190+
"ue3 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n",
191+
"ue = Tuple(ue1, ue2, ue3)\n",
192+
"\n",
193+
"u = element_of(V, name='u')\n",
194+
"\n",
195+
"error = Matrix([u[0]-ue[0], u[1]-ue[1], u[2]-ue[2]])\n",
196+
"\n",
197+
"# create the formal Norm object\n",
198+
"l2norm = Norm(error, domain, kind='l2')\n",
199+
"\n",
200+
"# discretize the norm\n",
201+
"l2norm_h = discretize(l2norm, domain_h, Vh)\n",
202+
"\n",
203+
"# assemble the norm\n",
204+
"l2_error = l2norm_h.assemble(u=uh)\n",
205+
"\n",
206+
"# print the result\n",
207+
"print(l2_error)"
208+
]
209+
},
210+
{
211+
"cell_type": "markdown",
212+
"id": "a6cbfeae",
213+
"metadata": {},
214+
"source": [
215+
"### Computing the $H^1$ semi-norm"
216+
]
217+
},
218+
{
219+
"cell_type": "code",
220+
"execution_count": null,
221+
"id": "e5c1a8b8",
222+
"metadata": {},
223+
"outputs": [],
224+
"source": [
225+
"# create the formal Norm object\n",
226+
"h1norm = SemiNorm(error, domain, kind='h1')\n",
227+
"\n",
228+
"# discretize the norm\n",
229+
"h1norm_h = discretize(h1norm, domain_h, Vh)\n",
230+
"\n",
231+
"# assemble the norm\n",
232+
"h1_error = h1norm_h.assemble(u=uh)\n",
233+
"\n",
234+
"# print the result\n",
235+
"print(h1_error)"
236+
]
237+
},
238+
{
239+
"cell_type": "markdown",
240+
"id": "3c09131c",
241+
"metadata": {},
242+
"source": [
243+
"### Computing the $H^1$ norm"
244+
]
245+
},
246+
{
247+
"cell_type": "code",
248+
"execution_count": null,
249+
"id": "d829e410",
250+
"metadata": {},
251+
"outputs": [],
252+
"source": [
253+
"# create the formal Norm object\n",
254+
"h1norm = Norm(error, domain, kind='h1')\n",
255+
"\n",
256+
"# discretize the norm\n",
257+
"h1norm_h = discretize(h1norm, domain_h, Vh)\n",
258+
"\n",
259+
"# assemble the norm\n",
260+
"h1_error = h1norm_h.assemble(u=uh)\n",
261+
"\n",
262+
"# print the result\n",
263+
"print(h1_error)"
264+
]
265+
}
266+
],
267+
"metadata": {
268+
"language_info": {
269+
"name": "python"
270+
}
271+
},
272+
"nbformat": 4,
273+
"nbformat_minor": 5
274+
}

0 commit comments

Comments
 (0)