diff --git a/examples/Freefem/MMS/1D/README.md b/examples/Freefem/MMS/1D/README.md index 84da35e..beb4ffa 100644 --- a/examples/Freefem/MMS/1D/README.md +++ b/examples/Freefem/MMS/1D/README.md @@ -54,6 +54,10 @@ The rule is determined by the choice of $n_g$, $\xi_k$, and $w_k$: |------|--------|----------|--------| | 1-point (midpoint) | 1 | $0$ | $2$ | | 2-point | 2 | $\pm \frac{1}{\sqrt{3}}$ | $1$ | +| 3-point | 3 | $0,\ \pm\sqrt{\frac{3}{5}}$ | $\frac{8}{9}\ (k=2),\quad \frac{5}{9}\ (k=1,3)$ | + +> **3-point rule:** $\xi_1 = -\sqrt{\tfrac{3}{5}},\ \xi_2 = 0,\ \xi_3 = \sqrt{\tfrac{3}{5}}$; +> $w_1 = w_3 = \tfrac{5}{9},\ w_2 = \tfrac{8}{9}$ ### Body Force @@ -206,3 +210,48 @@ With a 2-point quadrature rule: $$ F_0 = \frac{h}{6}\left[2f(x_0) + f(x_1)\right], \quad F_i = h\, f(x_i), \quad F_N = \frac{h}{6}\left[f(x_{N-1}) + 2f(x_N)\right] $$ + +$$ +F_i^{(e)} = \frac{h}{2} \sum_{k=1}^{2} w_k\, f(x_k^{(e)})\, \phi_i(x_k^{(e)}) +$$ + +--- +## Case 4: Exponential Function + +### Manufactured Solution + + +$$u_{ex}(x) = e^x - 1$$ + +### Source Term + +$$ +f(x)=−E exf(x) = -E\, e^xf(x)=−Eex +$$ + +### Boundary Conditions + + +$$ +u(0)=0 +$$ + + +$$ +\left.\frac{du}{dx}\right|_{1} = e +$$ + +Neumann force: + +$$ +F_N = E\, e +$$ + + +### Source Term Discretization + +With a 3-point quadrature rule: + +$$ +F_i^{(e)} = \frac{h}{2} \sum_{k=1}^{3} w_k\, f(x_k^{(e)})\, \phi_i(x_k^{(e)}) +$$ \ No newline at end of file diff --git a/examples/Freefem/MMS/1D/exponential.py b/examples/Freefem/MMS/1D/exponential.py new file mode 100644 index 0000000..1d44a26 --- /dev/null +++ b/examples/Freefem/MMS/1D/exponential.py @@ -0,0 +1,38 @@ +""" +Exponential MMS (non-dimensional, x = x_dim/L ∈ [0,1], E = E_dim/L): + + u_ex(x) = exp(x) - 1 du/dx = exp(x) + + f(x) = -E * exp(x) +BC: + u(0) = 0 (Dirichlet) + u'(1) = e (Neumann) => F_N = E·u'(1) = E·e +""" +from manufactured_solution import MMSCase1D +from bar import load_params, line_quadrature, build_bar_scene, run_scene +import numpy as np + +class Exponential(MMSCase1D): + name = "exponential" + plot_label = r"$e^x - 1$" + quadrature = staticmethod(line_quadrature(3)) # using 3 pts quadrature + + def u_ex(self, xi): + return np.exp(xi) - 1.0 + + def du_ex(self, xi): + return np.exp(xi) + + def source(self, xi, E): + return -E * np.exp(xi) # f = -E·u'' = -E·exp(x) + +mms = Exponential() + +def createScene(rootNode): + cfg = load_params() + L, E = cfg["length"], cfg["youngModulus"] + build_bar_scene(rootNode, mms, E / L, cfg["nx"]) + return rootNode + +if __name__ == "__main__": + run_scene(mms) \ No newline at end of file diff --git a/examples/Freefem/MMS/1D/params.json b/examples/Freefem/MMS/1D/params.json index 3aaf25e..4f354c2 100644 --- a/examples/Freefem/MMS/1D/params.json +++ b/examples/Freefem/MMS/1D/params.json @@ -4,7 +4,8 @@ "nx": 10, "nxConvergence": { "quadratic": [2, 3, 4, 5, 6, 8, 10, 12, 16, 20, 32, 64], - "cubic": [2, 3, 4, 5, 6, 8, 10, 12, 16, 20, 32, 64], - "sinusoidal": [10, 20, 40, 80, 160] + "cubic": [3, 4, 5, 6, 8, 10, 12, 16, 20, 32, 64], + "sinusoidal": [10, 20, 40, 80, 160], + "exponential" :[10, 20, 40, 80, 160] } } diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py index a00cd37..f7f58a4 100644 --- a/examples/Freefem/MMS/1D/run_convergence.py +++ b/examples/Freefem/MMS/1D/run_convergence.py @@ -7,6 +7,7 @@ from quadratic import mms as quadratic_mms from cubic import mms as cubic_mms from sinusoidal import mms as sinusoidal_mms +from exponential import mms as exponential_mms from bar import ( RESULTS_DIR, @@ -112,5 +113,5 @@ def run(mms, L, E, nx_list): if __name__ == "__main__": cfg = load_params() L, E = cfg["length"], cfg["youngModulus"] - for mms in (quadratic_mms, cubic_mms, sinusoidal_mms): + for mms in (quadratic_mms, cubic_mms, sinusoidal_mms, exponential_mms): run(mms, L, E, cfg["nxConvergence"][mms.name]) diff --git a/examples/Freefem/MMS/fem.py b/examples/Freefem/MMS/fem.py index 3bd692b..b73feb6 100644 --- a/examples/Freefem/MMS/fem.py +++ b/examples/Freefem/MMS/fem.py @@ -15,6 +15,8 @@ np.array([2.0])), 2: (np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)]), np.array([1.0, 1.0])), + 3: (np.array([-np.sqrt(3.0 / 5.0), 0.0, np.sqrt(3.0 / 5.0)]), + np.array([5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0])), }