<jupyter_start><jupyter_text>Introduction to IGA Finite Elements methodLet $\Omega \subset \mathbb{R}^d$ be a computational domain that is the image of a logical domain $\mathcal{P}$, *i.e.* a unit line (in *1d*), square (in *2d*) or a cube (in *3d*) with a **mapping** function $$\Omega = \mathcal{G} (\mathcal{P}) $$We consider the Poisson problem with Homogeneous Dirichlet boundary conditions:$$  \mbox{Find} ~  u \in H^1_0(\Omega), ~ \mbox{such that}$$$$  - \nabla^2 u = f, ~~ \Omega$$Using a Galerkin-Rietz (or more generally Galerkin-Petrov) method, we introduce:1. a discrete finite elements space $\mathcal{V}_h = \mathbf{span}\{ \phi_j, 1 \leq j \leq n_V \}$ for trial functions2. a discrete finite elements space $\mathcal{W}_h = \mathbf{span}\{ \psi_i, 1 \leq i \leq n_W \}$ for test functions (in the case of Galerkin-Rietz method, we have  $\mathcal{W}_h = \mathcal{V}_h$ )---**Note**:For the Poisson problem, we only need the Rietz-Galerkin method. ------**Note**:Here the index $i$ is a multi-index. For example, in $2D$, we have $i = (i_1, i_2)$------**Note**:The basis functions for both the trial and test spaces, must be of compact support, and behave like polynomials**TODO cite <NAME>**---Let $u_h \in \mathcal{V}_h$ such that $u_h = \sum_{j=1}^{n_V} u_j \phi_j$. Then the weak formulation for the Poisson equation writes$$  \sum_{j=1}^{n_V} u_j \int_{\Omega} \nabla \phi_j \cdot \nabla \psi_i = \int_{\Omega} f \psi_i, \quad \forall  1 \leq i \leq n_W $$Now, because we are using the Galerkin-Rietz approach, we have$$  \sum_{j=1}^{n} u_j \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i = \int_{\Omega} f \phi_i, \quad \forall  1 \leq i \leq n $$which can be written in a matrix form$$  M U = F$$where $U$ denotes the coefficients $(u_j, ~ 1 \leq j \leq n)$ and $F$ is a vector consisting of the terms $\int_{\Omega} f \phi_i$ for $1 \leq i \leq n$. Finally, the matrix $M$ is called the **stiffness** matrix and its entries are defined as$$  M_{ij} = \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i$$We will denote our basis function by $b_i$ and $b_j$ rather than $\phi_i$ and $\phi_j$. In this case, in order to solve the Poisson problem, one needs to1. Assemble the matrix $M$ of entries$$  M_{ij} = \int_{\Omega} \nabla b_j \cdot \nabla b_i$$2. Assemble the right hand side $F$ of entries$$  F_{i} = \int_{\Omega} \nabla f b_i$$3. Solve the linear system $$  M U = F$$  Assembly procedureLet's take a subdivion of the unit domain $\mathcal{P}$ as a structured grid. We can write $\mathcal{P} = \cup_e Q_e $ where $Q_e \cap Q_f = \emptyset$ when $e \ne f$. This defines also a partition over the computational domain $\Omega = \cup_e K_e$ where $K_e = \mathcal{G} (Q_e)$, ---**Note**:For the moment, we will assume that the mapping is the identity function, which means that $Q_e = K_e$. We will get back to the general case later.------**Note**:As in the basis functions, the index $e$ denote a multi-index.---Now let's go back to a matrix entry $M_{ij}$, we have$$M_{ij} = \sum_e \int_{Q_e} \nabla b_i \cdot \nabla b_j$$We know that every basis function $b_i$ is a polynomial over the element $Q_i$. We can then use one of the Gauss quadrature formulae. ---**Note**:In the case of *B-splines* we will avoid having to evaluate on the boundary of our element (B-Splines are only defined on the right, if we don't do so, evaluating derivatives may lead to wrong results). Hence we will use the Gauss-Legendre formula.------**Note**:In general, the quadrature formulae are given on the interval $[-1, 1]$. We therefor need to remap these points for every element, in out partition.---Let's assume now that our quadrature points have been remapped, and let's call them $\xi_k$, $\omega_k$ will denote the associated weight to $\xi_k$.In this case $M_{ij}$ can be written as$$M_{ij} = \sum_e \sum_k \omega_k \nabla b_i(\xi_k) \cdot \nabla b_j( \xi_k)$$  Naive procedureThe first simplest approach, is to evaluate every basis function $b_i$ over all the grid and then do the double sum. This may work in $1D$, but will be difficult to extend to higher dimensions, since it will consume too much memory. In fact, we do not take advantage of the compact support property of the basis functions!---**Exercise**:Using the *bsplines.py* file, implement the naive assembly procedure in $1D$---  Local procedureSince we are considering a *B-splines* discretization, we should rewrite our assembly procedure and use some of their properties.- For every element $e = (e_1, \ldots, e_d)$, we will only evaluate the non-vanishing *B-splines*- The evaluation will be done only in $1D$ - We will do this evaluation as a *pre-process* and prepare the *B-splines* values (and their derivaties), gathered for every *1D* element.At the end, the assembly procedure will look like:<jupyter_code>def assemble_stiffness(ne1,n1,p1,k1,spans_1,basis_1,weights_1):
    """Assembles the stiffness matrix in 1D.
    
    ne1: int
        number of elements
    n1: int
        total number of bsplines
    p1: int
        spline degree
    k1: int
        number of quadrature points per element
    spans_1: array
        array given the index of the last non-vanishing spline 
        for every element
    basis_1: array
        values of non-vanishing B-splines and their derivatives 
        on the quadrature points for every elemetn
    weights_1: array
        associated weights for the quadrature points for every element,
        where the size of the current element is taken into account
    """
    # ... 
    M = zeros(n1, n1)
    # ...

    # ... build matrices
    for ie1 in range(0, ne1):
        i_span_1 = spans_1[ie1]
        for il_1 in range(0, p1+1):
            for jl_1 in range(0, p1+1):
                i1 = i_span_1 - p1  - 1 + il_1
                j1 = i_span_1 - p1  - 1 + jl_1

                v_s = 0.0
                for g1 in range(0, k1):
                    bi_0 = basis_1[il_1, 0, g1, ie1]
                    bi_x = basis_1[il_1, 1, g1, ie1]

                    bj_0 = basis_1[jl_1, 0, g1, ie1]
                    bj_x = basis_1[jl_1, 1, g1, ie1]

                    wvol = weights_1[g1, ie1]

                    v_s += (bi_x * bj_x) * wvol

                M[i1, j1]  += v_s
    # ...
    
    return M<jupyter_output><empty_output><jupyter_text>In fact, we are not done yet!The problem now, is that the data structure we are using for the matrix (here as a dense matrix) will consume too much memory, and is not taken into account the locality!Since on each element, there are exactly $p+1$ consecutive non-vanishing *B-splines*, we know that on each element there are at most $2p+1$ non zeros entries (in every direction). This means that we will need to set entries as the following, in $1D$:```pythonM[i1, j1 - i1]  += v_s```in $2D$:```pythonM[i1, i2, j1 - i1, j2 - i2]  += v_s```and in $3D$:```pythonM[i1, i2, i3, j1 - i1, j2 - i2, j3 - i3]  += v_s```Unfortunatly, this can not be represented as a **numpy.ndarray**, because of the negative indexing!Moreover, if our aim is to write a parallel code, other consideration must be taken into account.More details can be found in the next session.**TODO reference to the notebook**<jupyter_code># css style
from IPython.core.display import HTML
def css_styling():
    styles = open("../../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11910569105691057, 'fraction_numerical': 0.017208672086720868, 'mean_word_length': 3.181869688385269, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '321146', 'license': 'MIT', 'n_tokens_mistral': 2432, 'n_tokens_neox': 2343, 'n_words': 1019, 'path': 'lessons/Chapter0/02_introduction_fem.ipynb', 'repo_name': 'ratnania/MHDPython', 'size': '7384'}