.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_poisson.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_poisson.py: .. currentmodule:: fdg N-D Poisson =========== This example demonstrates system for N-dimensional mixed Poisson equation can be set up. Mixed Poisson equation is defined in the weak form as: .. math:: :label: examples_nd_poisson_1 \left( p^{(n - 1)}, q^{(n - 1)} \right)_\Omega + \left( \mathrm{d} p^{(n - 1)}, u^{(n)} \right)_\Omega = \int_{\partial \Omega} p^{(n - 1)} \wedge \star u^{(n)} .. math:: :label: examples_nd_poisson_2 \left( v^{(n)}, \mathrm{d} q^{(n - 1)} \right)_\Omega = \left( v^{(n)}, f^{(n)} \right)_\Omega .. GENERATED FROM PYTHON SOURCE LINES 24-50 .. code-block:: Python from time import perf_counter from typing import Sequence import matplotlib.pyplot as plt import numpy as np import numpy.typing as npt from fdg import ( BasisSpecs, BasisType, CoordinateMap, DegreesOfFreedom, FunctionSpace, IntegrationMethod, IntegrationSpace, IntegrationSpecs, KForm, KFormSpecs, SpaceMap, compute_kform_mass_matrix, incidence_kform_operator, projection_kform_l2_dual, reconstruct, transform_kform_to_target, ) .. GENERATED FROM PYTHON SOURCE LINES 51-71 The manufactured solution for the general N-dimensional case uses the following for the manufactured solution: .. math:: :label: examples_nd_poisson_man_sol u^{(n)}(x_1, \dots, x_n) = k \left(\prod\limits_{i=1}^n \cos\left( \frac{\pi}{2} x_i \right)\right) \mathrm{d} x_1 \wedge \dots \wedge \mathrm{d} x_n This gives the forcing function: .. math:: :label: examples_nd_poisson_man_for f^{(n)}(x_1, \dots, x_n) = - k n \left(\frac{\pi}{2}\right)^2 \left( \prod\limits_{i=1}^n \cos\left( \frac{\pi}{2} x_i \right)\right) \mathrm{d} x_1 \wedge \dots \wedge \mathrm{d} x_n .. GENERATED FROM PYTHON SOURCE LINES 72-93 .. code-block:: Python SCALE = 0.1 def manufactured_solution(*x: npt.NDArray[np.double]) -> npt.NDArray[np.double]: """Exact manufactured solution.""" res = np.cos(x[0] * np.pi / 2) for v in x[1:]: res *= np.cos(v * np.pi / 2) return res * SCALE def manufactured_source_poisson(*x: npt.NDArray[np.double]) -> npt.NDArray[np.double]: """Exact manufactured source term.""" res = np.cos(x[0] * np.pi / 2) for v in x[1:]: res *= np.cos(v * np.pi / 2) res *= -((np.pi / 2) ** 2) * len(x) return res * SCALE .. GENERATED FROM PYTHON SOURCE LINES 94-106 First the geometry of the space this will be solved on will be defined. For this case, we use the unit square, where the interior is deformed, while boundaries are the same. The mapping for each coordinate is based on Equation :eq:`examples_nd_poisson_deformation`, with :math:`c` being the parameter that determines the scale of deformation. .. math:: :label: examples_nd_poisson_deformation x_i = \xi_i + c \prod\limits_{j=1}^n \left( 1 - {x_j}^2 \right) \sin \pi x_j .. GENERATED FROM PYTHON SOURCE LINES 107-137 .. code-block:: Python def disturbed_mapping( c: float, idx: int, *x: npt.NDArray[np.double] ) -> npt.NDArray[np.double]: """Return a perturbed map, where the boundaries are not affected, but the interior is. Parameters ---------- c : float Strenght of the disturbance. idx : int Index of the input to base the mapping on. *x : array Coordinates where the mapping should be computed. Returns ------- array Input coordinate ``idx``, but somewhat. """ base = x[idx] d = np.full_like(base, c) for v in x: d *= (1 - v**2) * np.sin(np.pi * v) return base + d .. GENERATED FROM PYTHON SOURCE LINES 138-141 Mappings for every coordinate are collected into a joined :class:`SpaceMap`, which is then used to map :math:`k`-form components between the reference domain and the physical domain. .. GENERATED FROM PYTHON SOURCE LINES 142-166 .. code-block:: Python def create_space_map( c: float, orders: Sequence[int], space: IntegrationSpace ) -> SpaceMap: """Create space map that are is disturbed.""" func_space = FunctionSpace( *(BasisSpecs(BasisType.LAGRANGE_UNIFORM, order) for order in orders) ) ndim = len(orders) points = np.meshgrid( *[np.linspace(-1, +1, order + 1) for order in orders], indexing="ij" ) return SpaceMap( *[ CoordinateMap( DegreesOfFreedom(func_space, disturbed_mapping(c, idim, *points)), space, ) for idim in range(ndim) ] ) .. GENERATED FROM PYTHON SOURCE LINES 167-176 With these two utilities, we can start working on computing convergence of the FEM discretization in the :math:`L^2`-norm First we must define how integration will be done. Here two :class:`IntegrationSpace` objects are defined - one for computing our results and another, finer, for computing the error. With these discretizations defined, we can use previously written functions to create the space mappings. .. GENERATED FROM PYTHON SOURCE LINES 177-192 .. code-block:: Python def create_space_maps(order_integration, type_integration, ndim, dp): """Create integration spaces.""" int_space = IntegrationSpace( *((IntegrationSpecs(order_integration, type_integration),) * ndim) ) int_space_higher = IntegrationSpace( *((IntegrationSpecs(order_integration + dp, type_integration),) * ndim) ) space_map = create_space_map(0.1, [5] * ndim, int_space) space_map_high = create_space_map(0.1, [5] * ndim, int_space_higher) return space_map, space_map_high .. GENERATED FROM PYTHON SOURCE LINES 193-200 Along with the :class:`IntegrationSpace` objects to define integration we must define the discretization of the :math:`k`-forms using a :class:`FunctionSpace` object. With base function space defined, we can define some :math:`k`-form specifications using :class:`KFormSpecs`. These do not contain any degrees of freedom themselves, but provide information about the order and function spaces. .. GENERATED FROM PYTHON SOURCE LINES 201-212 .. code-block:: Python def create_kform_specs(type_basis, order_basis, ndim): """Create k-form specifications.""" base_space = FunctionSpace(*((BasisSpecs(type_basis, order_basis),) * ndim)) specs_u = KFormSpecs(ndim, base_space) specs_q = KFormSpecs(ndim - 1, base_space) return specs_u, specs_q .. GENERATED FROM PYTHON SOURCE LINES 213-221 While for left side symmetry could be exploited, there's no harm to compute it in full for this example. To that end, we first compute the two mass matrices for the two :math:`k`-forms, then apply the incidence operators as needed. Now that we have the required blocks, we can assemble the system matrix. Alternatively we could have used Schur's complement to compute the solution, but for the sake of simplicity here we use the full dense solve. .. GENERATED FROM PYTHON SOURCE LINES 221-244 .. code-block:: Python def assemble_lhs(sm, specs_q, specs_u): """Crate the system matrix.""" mq = compute_kform_mass_matrix( sm, specs_q.order, specs_q.base_space, specs_q.base_space ) mu = compute_kform_mass_matrix( sm, specs_u.order, specs_u.base_space, specs_u.base_space ) mu_e = incidence_kform_operator(specs_q, mu, right=True) et_mu = incidence_kform_operator(specs_q, mu, transpose=True) system_matrix = np.block( [ [mq, et_mu], [mu_e, np.zeros_like(mu)], ] ) return system_matrix .. GENERATED FROM PYTHON SOURCE LINES 245-247 The right side of the Poisson equation can be computed from the "dual projection" of the manufactured source term on the function space. .. GENERATED FROM PYTHON SOURCE LINES 247-261 .. code-block:: Python def assemble_rhs(specs_u, specs_q, sm_h): """Assemble the system's RHS.""" source_vals = projection_kform_l2_dual([manufactured_source_poisson], specs_u, sm_h)[ 0 ] rhs = np.concatenate( (np.zeros(sum(specs_q.component_dof_counts)), source_vals.flatten()) ) return rhs .. GENERATED FROM PYTHON SOURCE LINES 262-266 From here we can split solution vector into degrees of freedom of individual :math:`k`-forms represented by :class:`KForm` objects. To compute the :math:`L^2` error, we need to reconstruct the computed solution, subtract the manufactured solution, then integrate the square of the error. .. GENERATED FROM PYTHON SOURCE LINES 266-298 .. code-block:: Python def reconstruct_error_l2(specs_q, specs_u, solution_dofs, ndim, sm_h): """Reconstruct the solution and compute the L2 error.""" sol_q = KForm(specs_q) sol_u = KForm(specs_u) nq = sum(specs_q.component_dof_counts) sol_q.values[:] = solution_dofs[:nq] sol_u.values[:] = solution_dofs[nq:] u_dofs = DegreesOfFreedom( specs_u.get_component_function_space(0), sol_u.get_component_dofs(0) ) # K-form computed values at integration nodes computed_values = transform_kform_to_target( ndim, sm_h, [reconstruct(u_dofs, *sm_h.integration_space.nodes())], )[0] # K-form exact values at integration nodes real_values = manufactured_solution( *[sm_h.coordinate_map(idx).values for idx in range(ndim)] ) # Error err_l2 = np.sum( (computed_values - real_values) ** 2 * sm_h.determinant * sm_h.integration_space.weights() ) return err_l2 .. GENERATED FROM PYTHON SOURCE LINES 299-301 All the small building blocks discussed before can now be put together to form the error calculation function. .. GENERATED FROM PYTHON SOURCE LINES 301-328 .. code-block:: Python def compute_l2_error( order_integration: int, type_integration: IntegrationMethod, order_basis: int, type_basis: BasisType, ndim: int, dp: int, ) -> float: """Solve the N-dimensional Poisson equation and compute the L^2 error.""" # Space maps sm, sm_h = create_space_maps(order_integration, type_integration, ndim, dp) # K-form specs specs_u, specs_q = create_kform_specs(type_basis, order_basis, ndim) # LHS of the system lhs = assemble_lhs(sm, specs_q, specs_u) # RHS rhs = assemble_rhs(specs_u, specs_q, sm_h) # Solve solution_dofs = np.linalg.solve(lhs, rhs) # Compute error squared err_l2 = reconstruct_error_l2(specs_q, specs_u, solution_dofs, ndim, sm_h) # Retur the error return float(np.sqrt(err_l2)) .. GENERATED FROM PYTHON SOURCE LINES 329-331 For this test, we will use Bernstein basis, Gauss integration rule, and the order difference of 1 between the lower and higher order integration rules. .. GENERATED FROM PYTHON SOURCE LINES 331-371 .. code-block:: Python BTYPE = BasisType.BERNSTEIN ITYPE = IntegrationMethod.GAUSS DP = 1 pvals = np.arange(1, 7) evals = np.zeros(pvals.size) tvals = np.zeros(pvals.size) for ndim in range(1, 4): for ip, p in enumerate(pvals): pv = int(p) t0 = perf_counter() l2 = compute_l2_error(pv + 0, ITYPE, pv, BTYPE, ndim, DP) t1 = perf_counter() evals[ip] = l2 tvals[ip] = t1 - t0 k1, k0 = np.polyfit(pvals, np.log(evals), deg=1) c = np.exp(k0) b = np.exp(k1) fig, ax = plt.subplots() ax.scatter(pvals, evals) ax.plot( pvals, c * b**pvals, linestyle="dashed", label=f"$\\varepsilon = {c:.2g} \\cdot {b:.2g}^{{p}}$", ) ax.set( yscale="log", xlabel="$p$", ylabel="$\\left|\\left| \\varepsilon \\right|\\right|_{ L^2 }$", ) ax.grid() ax.legend() ax.set_title(f"{ndim}-dimensional Poisson equation convergence") fig.tight_layout() plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/images/sphx_glr_plot_poisson_001.png :alt: 1-dimensional Poisson equation convergence :srcset: /auto_examples/images/sphx_glr_plot_poisson_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_poisson_002.png :alt: 2-dimensional Poisson equation convergence :srcset: /auto_examples/images/sphx_glr_plot_poisson_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/images/sphx_glr_plot_poisson_003.png :alt: 3-dimensional Poisson equation convergence :srcset: /auto_examples/images/sphx_glr_plot_poisson_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 8.235 seconds) .. _sphx_glr_download_auto_examples_plot_poisson.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_poisson.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_poisson.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_poisson.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_