Método dos Elementos Finitos

2 Problemas bidimensionais 2.2 Problema modelo 2.4 Malha adaptativa dual

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

2.3 Malha adaptativa baseada no resíduo

Em construção

Vamos estudar um método de malha adaptativa baseado no resíduo. Para isso, vamos precisar de uma estimativa a posteriori do erro.

2.3.1 Estimativa a posteriori

Teorema 2.3.1.(Estimativa a posteriori)

A solução $u_{h}$ do problema de elementos finitos (2.42) satisfaz

|||u-u_{h}|||^{2}\leq C\sum_{K\in\mathcal{K}}\eta_{K}^{2}(u_{h}),

(2.112)

onde o elemento residual $\eta_{K}(u_{h})$ é definido por

\eta_{K}(u_{h})=h_{K}\|f+\Delta u_{h}\|_{L^{2}(K)}+\frac{1}{2}h_{K}^{1/2}\|[% \boldsymbol{n}\cdot\nabla u_{h}]\|_{L^{2}(\partial K\setminus\partial\Omega)}.

(2.113)

Aqui, $[\boldsymbol{n}\cdot\nabla u_{h}]|_{K}$ denota o salto na derivada normal de $u_{h}$ nos lados interiores dos elementos de $\mathcal{K}$ . Além disso, lembremos que $\Delta u_{h}=0$ .

Demonstração.

Denotando $e:=u-u_{h}$ o erro entre a solução do problema forte e a solução de elementos finitos, temos

	$\displaystyle\\|e\\|_{E}^{2}=\\|\nabla e\\|_{L^{2}(\Omega)}^{2}$		(2.114)
	$\displaystyle\text{}\quad=\int_{\Omega}\nabla e\cdot\nabla e\,dx$		(2.115)
	$\displaystyle\text{}\quad=\int_{\Omega}\nabla e\cdot\nabla(e-\pi e)\,dx.$		(2.116)

Nesta última equação, temos usado a ortogonalidade de Galerkin (Teorema 2.2.2). Daí, temos

	$\displaystyle\int_{\Omega}\nabla e\cdot\nabla(e-\pi e)\,dx=\sum_{K\in\mathcal{% K}}\int_{K}\nabla e\cdot\nabla(e-\pi e)\,dx$		(2.117)
	$\displaystyle\text{}\quad=\sum_{K\in\mathcal{K}}-\int_{K}\Delta e(e-\pi e)\,dx$
	$\displaystyle\text{}\qquad+\int_{\partial K}\boldsymbol{n}\cdot\nabla e(e-\pi e% )\,ds,$		(2.118)
	$\displaystyle\text{}\quad=\sum_{K\in\mathcal{K}}\int_{K}(f+\Delta u_{h})(e-\pi e% )\,dx$
	$\displaystyle\text{}\qquad+\int_{\partial K\setminus\partial\Omega}\boldsymbol% {n}\cdot\nabla e(e-\pi e)\,ds,$		(2.119)

uma vez que $-\Delta e|_{K}=f+\Delta u_{h}|_{K}$ e, ambos, $e$ e $\pi e$ se anulam em $\partial\Omega$ .

Para computarmos o segundo termo do lado direito da ultima equação, observamos que o erro em lado $E$ recebe contribuições dos dois elementos $K^{\pm}$ que compartilham $E$ . Com isso, temos

	$\displaystyle\int_{\partial K^{+}\cap\partial K^{-}}\boldsymbol{n}\cdot\nabla e% (e-\pi e)\,ds=\int_{E}(\boldsymbol{n}^{+}\cdot\nabla e^{+}(e^{+}-\pi e^{+})$
	$\displaystyle\text{}\qquad+\boldsymbol{n}^{-}\cdot\nabla e^{-}(e^{-}-\pi e^{-}% ))\,ds,$		(2.120)

onde utilizamos a notação $v^{\pm}=v|_{K^{\pm}}$ . Lembremos que o erro $e$ é contínuo e, portanto, $(e^{+}-\pi e^{+})|_{E}=(e^{-}-\pi e^{-})|_{E}$ . Ainda, $\nabla u$ é contínuo, logo $(\boldsymbol{n}^{+}\cdot\nabla u^{+}+\boldsymbol{n}^{-}\cdot\nabla u^{-})|_{E}=0$ . Entretanto, $\nabla u_{h}|_{E}$ não é geralmente contínuo, sendo apenas constante por partes. Assim sendo e denotando o salto $[\boldsymbol{n}\cdot\nabla u_{h}]:=(\boldsymbol{n}^{+}\cdot\nabla u_{h}^{+}+% \boldsymbol{n}^{-}\cdot\nabla u_{h}^{-})$ , temos

	$\displaystyle\int_{E}(\boldsymbol{n}^{+}\cdot\nabla e^{+}(e-\pi e)+\boldsymbol% {n}^{-}\cdot\nabla e^{-}(e-\pi e))\,ds$
	$\displaystyle\text{}\quad=-\int_{E}[\boldsymbol{n}\cdot\nabla u_{h}](e-\pi e)% \,ds.$		(2.121)

Com isso, temos

\displaystyle\sum_{K\in\mathcal{K}}\int_{\partial K\setminus\partial\Omega}% \boldsymbol{n}\cdot\nabla e(e-\pi e)\,ds=-\sum_{E\in\mathcal{E}_{I}}\int_{E}[n% \cdot\nabla u_{h}](e-\pi e)\,ds,

(2.122)

onde $\mathcal{E}_{I}$ é o conjunto dos lados interiores na triangularização $\mathcal{K}$ . Logo, retornando a (2.119), obtemos

	$\displaystyle\\|e\\|_{E}^{2}=\sum_{K\in\mathcal{K}}\int_{K}(f+\nabla u_{h})(e-% \pi e)\,dx$
	$\displaystyle\text{}\quad-\frac{1}{2}\int_{\partial K\setminus\partial\omega}[% n\cdot\nabla u_{h}](e-\pi e)\,ds.$		(2.123)

Nos resta, agora, estimarmos estes dois termos do lado direito.

A estimativa do primeiro, segue da desigualdade de Cauchy-Schwarz seguida da estimativa padrão do erro de interpolação, i.e.

	$\displaystyle\int_{K}(f+\Delta u_{h})(e-\pi e)\,dx\leq\\|f+\Delta u_{h}\\|_{L^{2% }(\Omega)}\\|e-\pi e\\|_{L^{2}(\Omega)}$		(2.124)
	$\displaystyle\text{}\qquad\leq\\|f+\Delta u_{h}\\|_{L^{2}(\Omega)}Ch_{K}\\|De\\|_{% L^{2}(\Omega)}$		(2.125)

Para estimarmos as contribuições dos lados, usamos a desigualdade do Traço [Larson2023a, Subseção 7.1.6]

\|v\|_{L^{2}(\partial K)}^{2}\leq C\left(h_{K}^{-1}\|v\|_{L^{2}(K)}^{2}+h_{K}% \|\nabla v\|_{L^{2}(K)}^{2}\right).

(2.126)

Com esta, a desigualdade de Cauchy-Schwarz e a estimativa padrão do erro de interpolação, temos

	$\displaystyle\int_{\partial K\setminus\partial\Omega}[n\cdot\nabla u_{h}](e-% \pi e)\,ds\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}\\|e-\pi e\\|_{L^{2}(% \partial K)}$		(2.127)
	$\displaystyle\text{}\quad\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}$
	$\displaystyle\text{}\qquad\cdot C\left(h_{K}^{-1}\\|e-\pi e\\|_{L^{2}(K)}^{2}+h_% {K}\\|D(e-\pi e)\\|_{L^{2}(K)}^{2}\right)^{1/2}$		(2.128)
	$\displaystyle\text{}\quad\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}Ch_{K% }^{1/2}\\|De\\|_{L^{2}(K)}.$		(2.129)

Daí, a estimativa segue das (2.125) e (2.129). ∎

2.3.2 Malha adaptativa

Em construção

⬇

1from mpi4py import MPI

2from petsc4py import PETSc

3from dolfinx import mesh

4from dolfinx import fem

6# malha

7domain = mesh.create_unit_square(MPI.COMM_WORLD,

8 nx = 4, ny = 4)

11# problema mef

12import ufl

13from dolfinx import default_scalar_type

14from dolfinx.fem.petsc import LinearProblem, assemble_vector, create_vector

16# loop over refinements

17for i in range(3):

19 V = fem.functionspace(domain, ('P', 1))

21 # condição de contorno

22 import numpy as np

23 uD = fem.Function(V)

24 uD.interpolate(lambda x: np.full(x.shape[1], 0.))

26 tdim = domain.topology.dim

27 fdim = tdim - 1

28 domain.topology.create_connectivity(fdim, tdim)

29 boundary_facets = mesh.exterior_facet_indices(domain.topology)

30 boundary_dofs = fem.locate_dofs_topological(V, fdim,

31 boundary_facets)

32 bc = fem.dirichletbc(uD, boundary_dofs)

34 u = ufl.TrialFunction(V)

35 v = ufl.TestFunction(V)

37 # termo de fonte

38 x = ufl.SpatialCoordinate(domain)

39 f = 2 * ufl.pi**2 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])

41 # forma bilinear

42 a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

44 # forma linear

45 L = f * v * ufl.dx

47 # problema linear

48 problem = LinearProblem(a, L, bcs=[bc],

49 petsc_options_prefix="ex_mef2d_modelo")

50 uh = problem.solve()

52 # gráfico

53 import pyvista

54 from dolfinx import plot

56 # solução numérica

57 topology, cell_types, geometry = plot.vtk_mesh(V)

58 grid_uh = pyvista.UnstructuredGrid(topology, cell_types, geometry)

59 grid_uh.point_data["uh"] = uh.x.array.real

60 grid_uh.set_active_scalars("uh")

62 # plot

63 p = pyvista.Plotter()

64 p.add_mesh(grid_uh, show_edges=True, scalars="uh", cmap="viridis", label="uh")

65 # p.add_mesh(contours, color="white", label="uex_interp")

66 # p.add_point_labels(np.array(pts), [f"{val:.2f}" for val in vals], font_size=10, text_color="white", shape_opacity=0.5, point_color="white", point_size=1)

67 p.add_axes()

68 p.add_legend()

69 p.view_xy()

70 p.show()

72 # compute the L2 norm of the residual at each cell

73 W0 = fem.functionspace(domain, ("DG", 0))

74 w0 = ufl.TestFunction(W0)

76 # Strong residual of -Delta(u) = f

77 r = f + ufl.div(ufl.grad(uh))

78 r2_form = fem.form((r**2) * w0 * ufl.dx)

80 r2_vec = create_vector(fem.extract_function_spaces(r2_form))

81 with r2_vec.localForm() as loc_r2:

82 loc_r2.set(0)

83 assemble_vector(r2_vec, r2_form)

84 r2_vec.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,

85 mode=PETSc.ScatterMode.REVERSE)

86 r2_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES,

87 mode=PETSc.ScatterMode.FORWARD)

89 # cell_residual_l2 = fem.Function(W0)

90 # cell_residual_l2.x.array[:] = np.sqrt(np.maximum(r2_vec.array, 0.0))

91 # cell_residual_l2.x.scatter_forward()

93 # print norm of cell residuals

94 # print(f"Refinement {i}: max cell residual L2 norm = {cell_residual_l2.x.array.max()}")

95 print(f"Refinement {i}: max cell residual L2 norm = {r2_vec.array.max()}")

97 # plot cell residuals (DG0 is cellwise, so use mesh topology + cell_data)

98 topology, cell_types, geometry = plot.vtk_mesh(domain, domain.topology.dim)

99 grid_residual = pyvista.UnstructuredGrid(topology, cell_types, geometry)

100 # grid_residual.cell_data["cell_residual_l2"] = cell_residual_l2.x.array.real

101 grid_residual.cell_data["cell_residual_l2"] = r2_vec.array.real

102 grid_residual.set_active_scalars("cell_residual_l2")

103

104 pl = pyvista.Plotter()

105 pl.add_mesh(grid_residual, show_edges=True, scalars="cell_residual_l2", cmap="inferno", label="cell_residual_l2")

106 pl.add_axes()

107 pl.add_legend()

108 pl.view_xy()

109 pl.show()

110

111 # mark cells for refinement based on residuals

112 # DG0 dof numbering is not guaranteed to match cell numbering.

113 tdim = domain.topology.dim

114 num_local_cells = domain.topology.index_map(tdim).size_local

115 local_cells = np.arange(num_local_cells, dtype=np.int32)

116 cell_eta = np.empty(num_local_cells, dtype=np.float64)

117 for c in local_cells:

118 dof = W0.dofmap.cell_dofs(c)[0]

119 # cell_eta[c] = cell_residual_l2.x.array[dof].real

120 cell_eta[c] = r2_vec.array[dof].real

121

122 local_max = cell_eta.max() if num_local_cells > 0 else 0.0

123 global_max = domain.comm.allreduce(local_max, op=MPI.MAX)

124 threshold = 0.5 * global_max

125 cells_to_refine = local_cells[cell_eta > threshold]

126 print(f"Refinement {i}: marking {len(cells_to_refine)} cells for refinement")

127

128 # # global mesh refinement

129 # from dolfinx.mesh import refine

130 # refine_out = refine(domain)

131 # domain = refine_out[0] if isinstance(refine_out, tuple) else refine_out

132

133 # local mesh refinement

134 from dolfinx.mesh import refine

135 if len(cells_to_refine) > 0:

136 marker_entities = cells_to_refine

137 try:

138 # Some dolfinx versions expect edges, not cells, for local refine.

139 marker_entities = mesh.compute_incident_entities(domain.topology,

140 cells_to_refine,

141 tdim, 1)

142 except Exception:

143 pass

144 refine_out = refine(domain, marker_entities)

145 domain = refine_out[0] if isinstance(refine_out, tuple) else refine_out

2.3.3 Exercícios

Em construção

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

	$\displaystyle\int_{\partial K\setminus\partial\Omega}[n\cdot\nabla u_{h}](e-% \pi e)\,ds\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}\\|e-\pi e\\|_{L^{2}(% \partial K)}$		(2.127)
	$\displaystyle\text{}\quad\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}$
	$\displaystyle\text{}\qquad\cdot C\left(h_{K}^{-1}\\|e-\pi e\\|_{L^{2}(K)}^{2}+h_% {K}\\|D(e-\pi e)\\|_{L^{2}(K)}^{2}\right)^{1/2}$		(2.128)
	$\displaystyle\text{}\quad\leq\\|[n\cdot\nabla u_{h}]\\|_{L^{2}(\partial K)}Ch_{K% }^{1/2}\\|De\\|_{L^{2}(K)}.$		(2.129)

Política de uso de dados

Política de uso de dados

Método dos Elementos Finitos

2.3 Malha adaptativa baseada no resíduo

2.3.1 Estimativa a posteriori

Teorema 2.3.1.(Estimativa a posteriori)

Demonstração.

2.3.2 Malha adaptativa

2.3.3 Exercícios

Envie seu comentário