""" In this script, we implement the hdMSEM for the Poisson problem with a manufactured solution in the crazy_mesh. ⭕ To access the source code, click on the [source] button at the right side or click on :download:`[Poisson_problem_hd.py]`. Dependence may exist. In case of error, check import and install required packages or download required scripts. © mathischeap.com """ from numpy import pi, sin, cos import numpy as np from crazy_mesh_hybrid import CrazyMeshHybrid, CrazyMeshHybridGlobalNumbering, CrazyMeshHybridLocalBoundaryDOFs from mimetic_basis_polynomials import MimeticBasisPolynomials from incidence_matrices import E_div from trace_matrices import TF from mass_matrices import MassMatrices from projection import Reduction from scipy import sparse as spspa from scipy.sparse import linalg as spspalinalg from assembly import assemble from L2_error import L2Error from L2_error_dual import L2ErrorDual # the manufactured solutions def phi_exact(x, y, z): return sin(2 * pi * x) * sin(2 * pi * y) * sin(2 * pi * z) def u_exact(x, y, z): return 2 * pi * cos(2 * pi * x) * sin(2 * pi * y) * sin(2 * pi * z) def v_exact(x, y, z): return 2 * pi * sin(2 * pi * x) * cos(2 * pi * y) * sin(2 * pi * z) def w_exact(x, y, z): return 2 * pi * sin(2 * pi * x) * sin(2 * pi * y) * cos(2 * pi * z) def f_exact(x,y,z): return 12 * pi**2 * sin(2 * pi * x) * sin(2 * pi * y) * sin(2 * pi * z) def zero(x, y, z): # div u + f = 0 return 0 * x * y * z def Poisson_hd(K, N, c, save=False): """ :param int K: We use a crazy mesh of :math:`K^3` elements. The domain decomposition is based on this crazy mesh. :param int N: We use mimetic polynomials of degree :math:`N`. :param float c: The deformation factor of the crazy mesh is :math:`c,\ 0\\leq c\\leq 0.25`. :param save: Bool. If we save the coefficients of the variables. :return: A tuple of several outputs: - The :math:`L^2\\text{-error}` of solution :math:`\\boldsymbol{u}^h`. - The :math:`H(\\mathrm{div})\\text{-error}` of solution :math:`\\boldsymbol{u}^h`. - The :math:`L^2\\text{-error}` of solution :math:`\\varphi^h`. - The :math:`L^2\\text{-error}` of the projection, :math:`f^h`. - The :math:`L^2\\text{-error}` of :math:`\\nabla\\cdot\\boldsymbol{u}^h+f^h`. - The :math:`L^\\infty\\text{-error}` of :math:`\\nabla\\cdot\\boldsymbol{u}^h+f^h`. - The :math:`\\widetilde{H}^1\\text{-error}` of solution :math:`\\varphi^h`. :example: >>> K = 2 >>> N = 3 >>> c = 0 >>> Poisson_hd(K, N, c) # doctest: +ELLIPSIS hdMSEM L^2-error of u^h: 0.1535... """ K = int(K) N = int(N) c1000 = int(c*1000) # define the crazy mesh ... crazy_mesh = CrazyMeshHybrid(c, K) # generate the global numbering (gathering matrix) and find boundary dofs. GM_crazy_mesh = CrazyMeshHybridGlobalNumbering(K, N) GM_TF = GM_crazy_mesh.TF BD_crazy_mesh = CrazyMeshHybridLocalBoundaryDOFs(K, N) B_dofs_FP_dict = BD_crazy_mesh.FP B_dofs_FP = dict() for bn in B_dofs_FP_dict: if bn != 'x_minus': dofs_on_side = B_dofs_FP_dict[bn] for m in dofs_on_side: if m not in B_dofs_FP: B_dofs_FP[m] = list() B_dofs_FP[m].extend(dofs_on_side[m]) B_dofs_TF_dict = BD_crazy_mesh.TF B_dofs_TF_EN = dict() B_dofs_TF_NA = dict() for bn in B_dofs_TF_dict: if bn == 'x_minus': dofs_on_side = B_dofs_TF_dict[bn] for m in dofs_on_side: if m not in B_dofs_TF_EN: B_dofs_TF_EN[m] = list() B_dofs_TF_EN[m].extend(dofs_on_side[m]) else: dofs_on_side = B_dofs_TF_dict[bn] for m in dofs_on_side: if m not in B_dofs_TF_NA: B_dofs_TF_NA[m] = list() B_dofs_TF_NA[m].extend(dofs_on_side[m]) # define the basis functions _bfN_ = 'Lobatto-' + str(N) mbf = MimeticBasisPolynomials(_bfN_, _bfN_, _bfN_) # generate incidence matrix, trace matrices and mass matrices E = E_div(N, N, N) T = TF(N, N, N) MF = list() for k in range(K): for j in range(K): for i in range(K): ct = crazy_mesh.CT_of_element_index(i, j, k) MM = MassMatrices(mbf, ct) MF.append(MM.FP) # reduction of source term f_exact, and u_exact f_exact_local = list() u_exact_local = list() f_L2 = list() for k in range(K): for j in range(K): for i in range(K): ct = crazy_mesh.CT_of_element_index(i, j, k) RD = Reduction(mbf, ct) f_dofs_local = RD.VP(f_exact) f_exact_local.append(f_dofs_local) L2e = L2Error(mbf, ct) f_L2_local = L2e.VP(f_dofs_local, f_exact) f_L2.append(f_L2_local**2) _temp_ = RD.FP((u_exact, v_exact, w_exact)) u_exact_local.append(spspa.csc_matrix( _temp_[:,np.newaxis])) f_L2 = np.sum(f_L2)**0.5 # generate local systems invA_List = list() g_List = list() B_List = list() RSA = list() # reduced systems Ax = b RSb = list() # reduced systems Ax = b num_basis_FP = 3*(N+1)*N**2 num_basis_VP = N**3 num_basis_TF = 6 * N**2 for k in range(K): for j in range(K): for i in range(K): m = i + j * K + k * K**2 A00 = MF[m] A01 = E.T A02 = - T.T A10 = E A11 = spspa.csc_matrix((num_basis_VP, num_basis_VP)) A12 = spspa.csc_matrix((num_basis_VP, num_basis_TF)) A20 = A02.T.tolil() A21 = spspa.csc_matrix((num_basis_TF, num_basis_VP)) A22 = spspa.lil_matrix((num_basis_TF, num_basis_TF)) b0 = spspa.csc_matrix((num_basis_FP, 1)) b1 = spspa.csc_matrix(- f_exact_local[m][:,np.newaxis]) b2 = spspa.lil_matrix((num_basis_TF, 1)) # now apply local boundary condition if m in B_dofs_TF_EN: EN_dofs = B_dofs_TF_EN[m] A20[EN_dofs, :] = 0 A22[EN_dofs, EN_dofs] = 1 if m in B_dofs_TF_NA: NA_dofs = B_dofs_TF_NA[m] A20[NA_dofs, :] = 0 A20[NA_dofs, B_dofs_FP[m]] = 1 b2[NA_dofs] = u_exact_local[m][B_dofs_FP[m]] A = spspa.bmat(([A00, A01], [A10, A11]), format='csc') invA = spspalinalg.inv(A) invA_List.append(invA) B = spspa.bmat(([A02],[A12]), format='csc') B_List.append(B) C = spspa.bmat(([A20, A21],), format='csc') D = A22.tocsc() g = spspa.bmat(([b0],[b1]), format='csc') g_List.append(g) h = b2.tocsc() __ = C @ invA rsA = D - __ @ B rsb = h - __ @ g RSA.append(rsA) RSb.append(rsb) # assemble the reduce systems to a global system A = assemble(RSA, GM_TF, GM_TF) b = assemble(RSb, GM_TF) # solve the global system using the direct solver provided by scipy lamb = spspalinalg.spsolve(A, b) # solve Ax=b, obtain x del A, b lambda_local = lamb[GM_TF] u_local = list() phi_local = list() for k in range(K): for j in range(K): for i in range(K): m = i + j * K + k * K**2 u_phi_local_m = invA_List[m] @ (g_List[m].toarray().ravel('F') - B_List[m] @ lambda_local[m]) u_local_m = u_phi_local_m[:num_basis_FP] phi_local_m = u_phi_local_m[num_basis_FP:] u_local.append(u_local_m) phi_local.append(phi_local_m) u_local = np.array(u_local) phi_local = np.array(phi_local) div_u_local = (- E @ u_local.T).T div_u_plus_f = (E @ u_local.T + np.array(f_exact_local).T).T dual_gradient = (-E.T, T.T) dG_phi = list() for k in range(K): for j in range(K): for i in range(K): m = i + j * K + k * K**2 dG_phi_m = dual_gradient[0] @ phi_local[m] + dual_gradient[1] @ lambda_local[m] dG_phi.append(dG_phi_m) dG_phi = np.array(dG_phi) u_L2 = list() div_u_L2 = list() phi_L2 = list() div_L2 = list() div_Linf = list() phi_dH1 = list() dH1_error = list() for k in range(K): for j in range(K): for i in range(K): m = i + j * K + k * K**2 ct = crazy_mesh.CT_of_element_index(i, j, k) L2E = L2Error(mbf, ct) L2Ed = L2ErrorDual(mbf, ct) u_l2_m = L2E.FP(u_local[m], (u_exact, v_exact, w_exact)) phi_dH1_m = L2Ed.FP(dG_phi[m], (u_exact, v_exact, w_exact)) dH1_error_m = L2Ed._FP_diff_(dG_phi[m], u_local[m], (zero, zero, zero)) div_u_l2_m = L2E.VP(div_u_local[m], f_exact) phi_l2_m = L2Ed.VP(phi_local[m], phi_exact) div_L2_m = L2E.VP(div_u_plus_f[m], zero) div_Linf_m = L2E.VP(div_u_plus_f[m], zero, n='infty') u_L2.append(u_l2_m**2) div_u_L2.append(div_u_l2_m**2) phi_L2.append(phi_l2_m**2) div_L2.append(div_L2_m**2) div_Linf.append(div_Linf_m) phi_dH1.append(phi_dH1_m**2) dH1_error.append(dH1_error_m**2) u_L2 = np.sum(u_L2)**0.5 div_u_L2 = np.sum(div_u_L2)**0.5 phi_L2 = np.sum(phi_L2)**0.5 div_L2 = np.sum(div_L2)**0.5 div_Linf = np.max(div_Linf) phi_dH1 = np.sum(phi_dH1)**0.5 dH1_error = np.sum(dH1_error)**0.5 u_Hdiv = np.sqrt(u_L2**2 + div_u_L2**2) phi_dH1 = np.sqrt(phi_L2**2 + phi_dH1**2) # print info and return print('hdMSEM') print("L^2-error of u^h: ", u_L2) print("H(div)-error of u^h: ", u_Hdiv) print("L^2-error of phi^h: ", phi_L2) print("L^2-error of projection f^h-f: ", f_L2) print("L^2-error of div(u^h)+f^h: ", div_L2) print("L^inf-error of div(u^h)+f^h: ", div_Linf) print("dual_H^1-error pf phi^h: ", phi_dH1) print("L^2-error of dual gradient phi^h minus u: ", dH1_error) if save: name_temp = f'results/hdMSEM_K{K}_N{N}_c{c1000}_' # noinspection PyTypeChecker np.savetxt(name_temp + 'u.txt', u_local) # noinspection PyTypeChecker np.savetxt(name_temp + 'phi.txt', phi_local) # noinspection PyTypeChecker np.savetxt(name_temp + 'du_plus_f.txt', div_u_plus_f) # noinspection PyTypeChecker np.savetxt(name_temp + 'lambda.txt', lambda_local) # noinspection PyTypeChecker np.savetxt(name_temp + 'dGphi.txt', dG_phi) return u_L2, u_Hdiv, phi_L2, f_L2, div_L2, div_Linf, phi_dH1, dH1_error if __name__ == '__main__': import doctest doctest.testmod() K = 2 N = 3 c = 0.25 Poisson_hd(K, N, c, save=False)