diff --git a/content/GA_2_2/README.md b/content/GA_2_2/README.md
index 600cceb6aa3636c14ce27f3b8c553eecae0dddfa..c6721630001a3a75862fc585d378077567fa350f 100644
--- a/content/GA_2_2/README.md
+++ b/content/GA_2_2/README.md
@@ -25,12 +25,13 @@ To run the code for this assignment **add a new Python package** to your `mude-b
conda activate mude-base
conda install ipympl
```
-Remember to install before running your notebook. If you did not do this, you may have to restart the kernel in your notebook.
+Remember to install _before_ running your notebook. If you did not do this, you may have to restart the kernel in your notebook.
-## Task Overview
+There are occasionally issues with the visualizations in VS Code. If you cannot see the animations, try uncommenting the cell magic lines (those beggining with `%matplotlib ...`). You can start a text search by using `CTRL+F`.
-...
+## Task Overview
+You may be able to divide work on the last two questions once the derivation and determination of boundary conditions is complete. Note that the first question in particular is essential exam practice and every group member should be able to complete it individually.
**End of file.**
diff --git a/content/GA_2_2/Report_solution.md b/content/GA_2_2/Report_solution.md
index efd5ad0373e2b33cf589410b868468dfb51a0874..484b6c1cef9a9021b4d6f466f5fb9e81cd669801 100644
--- a/content/GA_2_2/Report_solution.md
+++ b/content/GA_2_2/Report_solution.md
@@ -14,66 +14,300 @@ $$
or by including an image with your (clearly written) handwriting.
-### 1. Boundary conditions
+**Question 1: Derivation**
-*1a. What boundary conditions are actively enforced? On which part of the boundary are they enforced?*
+Follow the steps from strong form to discretized form to derive the expression $\mathbf{M}=\int_\Omega\mathbf{N}^T\mathbf{N}\,\mathrm{d}\Omega$ in the term $\mathbf{M}\dot{\mathbf{u}}$. You will only be assessed on how you deal with the term that contains the time derivative. The other terms exactly following the recipe outlined for the [Poisson equation in 2D](https://mude.citg.tudelft.nl/2024/book/fem/poisson2d.html) in the book.
-***(1)*** Dirichlet boundary conditions: $u=10$ on the bottom left edge
-*1b. On the remainder of the boundary, nothing is done in the implementation to enforce any boundary conditions. Give a mathematical expression for the boundary condition that is naturally applied. Also describe an observation about the obtained solution that confirms that this boundary condition is indeed satisfied.*
-***(1)*** Homogeneous Neumann boundary conditions: $\nabla u\cdot\mathbf{n}=0$.
-***(1)*** It can be observed that the gradient of the solution in perpendicular to the boundary is equal to zero. This is best observed from the 2D visualization, where iso-lines separating the different colours are normal to the boundary.
+**solution**
+Time dependent 2D Poisson equation:
-### 2. Integration scheme
-*2a. Which integration schemes are used to compute the element contributions to the $\mathbf{K}$ and $\mathbf{M}$ matrices? Comment on the locations and weights of the integration points.*
+$$
+\frac{\partial u}{\partial t} = \nu \left(\frac{\partial^{2} u}{{\partial x}^{2}} + \frac{\partial^{2} u}{{\partial y}^{2}}\right) + q
+$$
+
+$$
+\frac{\partial u}{\partial t} = \nu \nabla^{2} u + q
+$$
+
+1. Integrate and introduce weight functions
+
+$$
+\int_{\Omega} w \frac{\partial u}{\partial t} d\Omega = \int_{\Omega} w \nu \nabla^2 u d\Omega + \int_{\Omega} wq d\Omega
+$$
+
+2. Integration by parts: changing the second derivative, changing sign and introducing boundary condition (not necessary for the term with time-derivative)
+
+$$
+\int_{\Omega} w \frac{\partial u}{\partial t} d\Omega = - \int_{\Omega} \nu \nabla w \cdot \nabla u d\Omega + \int_{\Gamma} w \mu \nabla u \cdot \bar{n} d\Gamma + \int_{\Omega} w q d\Omega
+$$
+
+3. Substitute boundary condition (w= 0 on $\Gamma_D$ and $\mu \nabla u \cdot \mathbf{n}$ on $\Gamma_N$) (not necessary for the term with time-derivative):
+
+$$\int_{\Omega} w \frac{\partial u}{\partial t} d\Omega + \int_{\Omega} \nu \nabla w \cdot \nabla u d\Omega = \int_{\Gamma N} w h d\Gamma_N + \int_{\Omega} w q d\Omega$$
+
+4. Discretization:
+
+$$
+u^h = \mathbf{N}\mathbf{u}
+$$
+$$w^h = \mathbf{N}\mathbf{w}$$
+
+
+$$\nabla u^h = \mathbf{B}\mathbf{u}$$
+$$\nabla w^h = \mathbf{B}\mathbf{w}$$
+
+5. Substitute and take $\mathbf{u,w}$ out of the integral as they don't depend on $x$ and $y$
+
+$$
+\mathbf{w^T} \int_{\Omega} \mathbf{N}^T \mathbf{N} d\Omega \frac{\partial \mathbf{u}}{\partial t} + \mathbf{w^T} \int_{\Omega} \nu \mathbf{B}^T \mathbf{B} d\Omega \mathbf{u} = \mathbf{w^T} \int_{\Gamma \mathbf{N}} \mathbf{N}^T h d\Gamma_\mathbf{N} + \mathbf{w}^T \int_{\Omega} \mathbf{N}^T q d\Omega
+$$
+
+6. Eliminate $\mathbf{w}^T$
+
+$$ \int_{\Omega} \mathbf{N}^T \mathbf{N} d\Omega \frac{\partial \mathbf{u}}{\partial t}+ \int_{\Omega} \nu \mathbf{B}^T \mathbf{B} d\Omega \mathbf{u} = \int_{\Gamma \mathbf{N}} \mathbf{N}^T h d\Gamma_\mathbf{N} + \int_{\Omega} \mathbf{N}^T q d\Omega
+$$
+
+
+7. write as a system of equations:
+
+$$
+\mathbf{M} \frac{\partial \mathbf{u}}{\partial t} + \mathbf{K u} = \mathbf{q}$$
-***(1)*** For the $\mathbf{M}$-matrix: 3-points at the middle of the edges of the triangle, each with weight $A/3$ where $A$ is the element area
+$$ \mathbf{M \dot{u}} + \mathbf{K u}Â =Â \mathbf{q}$$
-***(1)*** For the $\mathbf{K}$-matrix there is no location because $\mathbf{K}$ does not depend on position, there is a single integration point with weight $A$ (or the integrand is just multiplied with $A$)
-
-*2b. Which changes would you make to the code to evaluate the $\mathbf{M}$-matrix with a single integration point at the center of gravity of each element (give your answer by indicating where you would replace existing code and typing out the new lines of code)*
-***(1)***
+**Question 2: Problem definition**
-In the function `get_element_M`, change the lines that define `integration_locations` and `integration_weights`:
+Investigate the code and results to find out which problem is being solved.
- `integration_locations = [(n1+n2+n3)/3]`
-
- `integration_weights = [element_area]`
+- Give a mathematical description of the problem in terms of governing equation and boundary conditions. Be as specific as possible, indicating the values that are used as input to the calculation.
+- In the final visualization contour lines are visible, connecting points that have the same temperature. As the solution evolves, these contour lines remain approximately perpendicular to the boundary. Which boundary condition does this observation relate to?
-### 3. Time step size dependence
-*3a. Try increasing the step size $\Delta t$. What is the reason this code does not suffer from instability for large time steps?*
+**solution**
-***(1)***
-The time-integration scheme is Euler backward.
+1. Governing equations and boundary conditions:
-*3b. Try decreasing the time step to very small numbers. If you make the time step small enough, some unphysical behavior can be observed in the solution, at least for initial time steps. What is the source of this behavior?*
-***(1)***
-For very small $t$, the solution initially has local high gradients. The mesh is not fine enough to resolve this accurately. With small $\Delta t$ the first time steps are in this domain.
+The governing equation is:
-### 4. $\mathbf{B}$-matrix
+$$
+\frac{\partial u}{\partial t} = \nu \Delta u + q,
+$$
+
+The input values are:
+
+- $\nu = 1$,
+- $q = 15$,
+
+Boundary conditions:
+
+Dirichlet boundary conditions: $u=10$ on the bottom right edge
+$$
+u(x, y, t) = 10, \quad \text{on } \Gamma_D,
+$$
+
+Homogeneous Neumann boundary conditions on remainder of the bondary
+$$
+\nabla u \cdot \mathbf{n} = 0, \quad \text{on } \Gamma_N,
+$$
+
+Initial conditions:
+$$
+u(x, y, 0) = 30, \quad \in \Omega.
+$$
+
+2. Which boundary condition does this observation relate to?
+
+- The homogeneous Neumman boundary condition
+
+**Question 3: Integration scheme**
+
+- In the `get_element_M` function, how many integration points are used and where in the triangles are they positioned?
-*4a Shape functions in the triangular element each have the form $N_i=a_ix+b_iy+c_i$ with $i\in[1,3]$. For every $i$, the coefficients $a_i, b_i, c_i$ are computed in the code to form the B-matrix. Why does the $\mathbf{B}$-matrix inside the element not depend on $x$ and $y$?*
+- In the `get_element_K` a simpler implementation is used. What is the essential difference between $\mathbf{K}_e$ and $\mathbf{M}_e$ that is the reason why this simpler implementation is valid for $\mathbf{K}_e$? (The subscript $_e$ is used to indicate the contribution to the matrix from a single element, or the *element matrix*).
+
+
+
+
+**solution**
+
+1.
+- 3 integration points
+- location at the midpoints of the triangle edges
+in the code:
+integration_locations = [(n1 + n2) / 2, (n2 + n3) / 2, (n3 + n1) / 2]
+
+2.
+
+$$ \mathbf{M} = \int_{\Omega} \mathbf{N}^T \mathbf{N} \,d \Omega$$
+
+$$ \mathbf{K} = \int_{\Omega} \mathbf{B}^T \nu \mathbf{B} \,d \Omega$$
+
+
+$$
+N_i = a_i + b_i x + c_i y, \quad i \in [1, 3],
+$$
+
+where $a_i, b_i, c_i$ are constants determined by the geometry of the triangle. The shape functions $N_i$ vary linearly across the element because they depend on $x$ and $y$.
+
+
+Taking the derivative of $N_i$ with respect to $x$ and $y$ gives:
+
+$$
+\frac{\partial N_i}{\partial x} = b_i, \quad \frac{\partial N_i}{\partial y} = c_i.
+$$
+
+The derivatives $b_i$ and $c_i$ are constants because $N_i$ is linear, and the derivative removes the dependence on $x$ and $y$.
-***(1)***
-The $\mathbf{B}$-matrix contains derivatives with respect to $x$ and $y$. Because there are only linear terms in $N_i$, these derivatives are not a function of $x$.
+The $\mathbf{B}$-matrix is therefore constant within a single element and does not vary with $x$ or $y$. This simplifies the computation of the stiffness matrix $\mathbf{K}$ because $\mathbf{B}^T \nu \mathbf{B}$ remains constant and only needs to be multiplied by the area of the triangle.
+
+
+**Question 4: Shape functions**
+Investigate the shape functions for the element with index 10 in the mesh. Use the `get_shape_functions_T3` function defined in the notebook to find expressions for the shape functions in that element and check that they satisfy the shape function properties.
+
+1. What are the coordinates of the nodes of element 10?
+
+2. What are the shape functions of the element?
+
+3. Assert that the shape functions satisfy the partition of unity property:
+
+$$
+\sum_i N_i(\mathbf{x}) = 1
+$$
-*4b Give an expression for the $\mathbf{B}$-matrix in terms of these nine coefficients ($a_1, a_2, a_3, b_1, b_2, b_3, c_1, c_2, c_3$).*
+4. Assert for one of the shape functions that it satisfies the Kronecker delta property
-***(1)***
$$
-\mathbf{B} = \left[\begin{matrix}
-a_1 & a_2 & a_3 \\
-b_1 & b_2 & b_3
-\end{matrix}\right]
+N_i(\mathbf{x}_j) = \begin{cases}
+ 1, & i = j \\
+ 0, & i\neq j
+\end{cases}
$$
+**solution**
+
+*1. Find coordinates of element 10*
+
+- `nodes, connectivity = readGmsh('bigM.msh')`
+
+`nodes` is an array of the coordinates of the nodes.
+
+`connectivity` is an array of the nodes belonging to each element.
+
+
+Node indices of element 10: [ 55, 56, 265]
+
+Coordinates of the nodes of element 10:
+| $N_i$ | x | y | z |
+|------------|--------------|--------------|--------------|
+| 1 | 0.1 | 1.0 | 0.0 |
+| 2 | 0.05 | 1.0 | 0.0 |
+| 3 | 0.07628505 | 0.94825032 | 0.0 |
+
+*2. shape functions*
+
+ The `get_shape_functions_T3` outputs the coefficients for each shape function in the form:
+
+ $$
+ N_i = a_i +b_i x+ c_iy
+ $$
+
+ These coeffients are stored in an 3x3 mtrix coeffs
+ These coeffients define the shape functions
+
+For the 3 nodes of element 10 we get the following coeffients:
+
+| Node | $a_i$ | $b_i$ | $c_i$ |
+|------------|--------------|--------------|--------------|
+| 1 | -11.15853885 | 20.0 | 10.15853885 |
+| 2 | -7.16525352 | -20.0 | 9.16525352 |
+| 3 | 19.32379237 | 0.0 | -19.32379237 |
+
+
+
+
+*3. Assert that the shape functions satisfy the partition of unity property*:
+
+$$
+\sum_{i=1}^3 N_i(\mathbf{x}) = 1
+$$
+
+ - Each shape function is of the form:
+ $$
+ N_i(x, y) = a_i + b_i x + c_i y, \quad i \in [1, 3],
+ $$
+ where $a_i, b_i, c_i$ are the coefficients computed using the `get_shape_functions_T3` function.
+
+ Choose a Point Inside the Element:
+ - for example centroid $\frac{N_1 +N_2+N_3}{3}$
+
+Evaluate the Shape Functions at the this Point.
+ - For each shape function $N_i$, compute its value at the test point
+
+ Compute the Sum of Shape Functions:
+ - Add the values of all shape functions at the test point:
+ $$\text{Sum} = \sum_{i=1}^3 N_i(\mathbf{x}).$$
+
+` sum(coeff[0] + coeff[1] * centroid[0] + coeff[2] * centroid[1] for coeff in coeffs)`
+
+
+ Assert the Partition of Unity:
+ - Check if the sum equals $1$:
+$$
+\text{Assert: } \sum_{i=1}^3 N_i(\mathbf{x}) = 1.
+$$
+
+
+```python
+# Get the coordinates of the nodes for element 10
+element_nodes = nodes[connectivity[9]]
+
+# Compute the shape function coefficients for element 10
+coeffs = get_shape_functions_T3(element_nodes[0], element_nodes[1], element_nodes[2])
+
+# Compute the centroid of the element
+centroid = element_nodes.mean(axis=0)
+
+# Evaluate the sum of the shape functions at the centroid
+partition_of_unity = sum(
+ coeff[0] + coeff[1] * centroid[0] + coeff[2] * centroid[1] for coeff in coeffs
+)
+
+# Print the result
+print(f"Sum of shape functions at the centroid: {partition_of_unity}")
+assert abs(partition_of_unity-1)< 1e-5,, "Partition of unity property not satisfied!"
+```
+
+*4. Assert for one of the shape functions that it satisfies the Kronecker delta property*
+
+For each shape function $N_i(x, y,z)$, calculate its value at the coordinates of each node $(x_j, y_j, z_j)$ then check if the shapefuction = 1
+
+```python
+# Extract coordinates of the nodes for element 10
+element_coords = nodes[connectivity[9]]
+
+# Compute shape function coefficients for element 10
+coeffs = get_shape_functions_T3(element_coords[0], element_coords[1], element_coords[2])
+
+# Verify Kronecker delta property
+for i in range(3):
+ for j in range(3):
+ x, y, z = element_coords[j]
+ N_i_at_xj = coeffs[i][0] + coeffs[i][1] * x + coefffs[i][2] * y
+ expected = 1.0 if i == j else 0.0
+ print(f"N_{i+1}(x_{j+1}) = {N_i_at_xj}, Expected = {expected}")
+ assert np.isclose(N_i_at_xj, expected), f"Kronecker delta property failed for N_{i+1} at node {j+1}"
+
+```
+
+
+
## General Comments on the Assignment [optional]
_Use this space to let us know if you encountered any issues completing this assignment (but please keep it short!). For example, if you encountered an error that could not be fixed in yout Python code, or perhaps there was a problem submitting something via GitLab. You can also let us know if the instructions were unclear. You can delete this section if you don't use it._
@@ -81,4 +315,4 @@ _Use this space to let us know if you encountered any issues completing this ass
**End of file.**
<span style="font-size: 75%">
-© Copyright 2024 <a rel="MUDE" href="http://mude.citg.tudelft.nl/">MUDE</a>, TU Delft. This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">CC BY 4.0 License</a>.
\ No newline at end of file
+© Copyright 2024 <a rel="MUDE" href="http://mude.citg.tudelft.nl/">MUDE</a>, TU Delft. This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">CC BY 4.0 License</a>.
diff --git a/content/GA_2_2/environment.yml b/content/GA_2_2/environment.yml
new file mode 100644
index 0000000000000000000000000000000000000000..765e2f0a7b1d3f844667c289566fb9c1d186e287
--- /dev/null
+++ b/content/GA_2_2/environment.yml
@@ -0,0 +1,11 @@
+name: mude-week-2-2
+dependencies:
+ - python=3.12
+ - numpy
+ - scipy
+ - pandas
+ - matplotlib
+ - jupyter
+ - ipykernel
+ - ipympl
+ - pip
\ No newline at end of file
diff --git a/content/GA_2_2/the_big_M.ipynb b/content/GA_2_2/the_big_M.ipynb
index f305f5c14f78c7973127d8fa3122987ba92d0fe2..f8540be29b5e03798eec41500081af4d4644acb6 100644
--- a/content/GA_2_2/the_big_M.ipynb
+++ b/content/GA_2_2/the_big_M.ipynb
@@ -170,7 +170,7 @@
" plt.axis('off')\n",
"\n",
"nodes, connectivity = readGmsh('bigM.msh')\n",
- "plotMesh(nodes, connectivity)"
+ "plotMesh(nodes, connectivity)\n"
]
},
{
@@ -402,7 +402,7 @@
"source": [
"# first option: plot a single time step in a plot that can be rotated manually\n",
"\n",
- "%matplotlib widget\n",
+ "# %matplotlib widget\n",
"from ipywidgets import interact, fixed, widgets\n",
"from matplotlib import colors\n",
"from matplotlib import cm\n",
@@ -447,7 +447,7 @@
"source": [
"# second option: make an animation to see the evolution\n",
"\n",
- "%matplotlib inline\n",
+ "# %matplotlib inline\n",
"\n",
"from ipywidgets import interact, fixed, widgets\n",
"from matplotlib import colors\n",
@@ -499,7 +499,7 @@
"source": [
"# third option: 2-dimensional representation\n",
"\n",
- "%matplotlib inline\n",
+ "# %matplotlib inline\n",
"\n",
"def plot_result(nodes, conn, results, step):\n",
" fig = plt.figure()\n",
@@ -557,19 +557,11 @@
"<span style=\"font-size: 75%\">\n",
"© Copyright 2024 <a rel=\"MUDE\" href=\"http://mude.citg.tudelft.nl/\">MUDE</a> TU Delft. This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">CC BY 4.0 License</a>."
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "694bc791",
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
"kernelspec": {
- "display_name": "mude-base",
+ "display_name": "mude-week-2-2",
"language": "python",
"name": "python3"
},
@@ -583,7 +575,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.12.5"
+ "version": "3.12.7"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,