Método dos Elementos Finitos

2 Problemas bidimensionais 2.3 Malha adaptativa baseada no resíduo 3 Aplicações 2D

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

2.4 Malha adaptativa dual

Em construção

Podemos desenvolver um método de malha adaptativa para minimizar o erro de uma quantidade de interesse, dada por um funcional linear $\mathcal{M}:V\to\mathbb{R}$ , onde $V$ é o espaço de soluções do problema fraco. Aqui, vamos seguir trabalhando com o problema modelo (2.37): encontrar $u\in V_{0}$ tal que

a(u,v)=L(v),\leavevmode\nobreak\ \forall v\in V_{0},

(2.130)

onde $V_{0}=H_{0}^{1}(\Omega)$ e $\Omega=[0,1]^{2}$ . Como exemplo de funcional linear, vamos usar a média da solução $u$ no domínio $\Omega$ , ou seja,

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@color@rgb@fill{0}{0}{1}\mathcal{M}(u)% =\frac{1}{|\Omega|}\int_{\Omega}u\,dx.

(2.131)

Tendo em vista que o funcional objetivo é linear, temos que

\mathcal{M}(e)=\mathcal{M}(u-u_{h})=\mathcal{M}(u)-\mathcal{M}(u_{h}),

(2.132)

onde $e=u-u_{h}$ é o erro entre a solução exata $u$ e a solução de elementos finitos $u_{h}$ . Introduzimos o problema dual (ou problema adjunto) associado a $\mathcal{M}$ : encontrar $z\in V_{0}$ tal que

a(v,z)=\mathcal{M}(v),\leavevmode\nobreak\ \forall v\in V_{0}.

(2.133)

Com isso, temos

	$\displaystyle\mathcal{M}(e)=a(e,z)$		(2.134)
	$\displaystyle\text{}\quad=a(e,z-z_{h})+a(e,z_{h})$		(2.135)
	$\displaystyle\text{}\quad=a(u-u_{h},z-z_{h})$		(2.136)
	$\displaystyle\text{}\quad=a(u,z-z_{h})-a(u_{h},z-z_{h})$		(2.137)
	$\displaystyle\text{}\quad=L(z-z_{h})-a(u_{h},z-z_{h}).$		(2.138)

Para o caso do problema modelo, temos

	$\displaystyle\mathcal{M}(e)=\int_{\Omega}f(z-z_{h})\,dx-\int_{\Omega}\nabla u_% {h}\cdot\nabla(z-z_{h})\,dx$		(2.139)
	$\displaystyle\text{}\quad=\sum_{K\in\mathcal{K}}\int_{K}(f+\Delta u_{h})(z-z_{% h})\,dx$
	$\displaystyle\text{}\qquad+\frac{1}{2}\sum_{E\in\mathcal{E}_{I}}\int_{E}[% \boldsymbol{n}\cdot\nabla u_{h}](z-z_{h})\,ds,$		(2.140)

que nos fornece uma forma de representação do erro.

Teorema 2.4.1.(Estimativa a posteriori dual)

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

\mathcal{M}(e)\leq C\sum_{K\in\mathcal{K}}\eta_{K}(u_{h},z),

(2.141)

onde o elemento residual $\eta_{K}(u_{h},z)=\rho_{K}(u_{h})\omega_{K}(z)$ é definido por

	$\displaystyle\rho_{K}(u_{h})=\\|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.142)
	$\displaystyle\omega_{K}(z)=h^{2}_{K}\\|D^{2}z\\|_{L^{2}(K)}.$		(2.143)

Código 17: adapt_dual.py

⬇

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(4):

19 # primal problem

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

22 ## condição de contorno

23 import numpy as np

24 uD = fem.Function(V)

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

27 tdim = domain.topology.dim

28 fdim = tdim - 1

29 domain.topology.create_connectivity(fdim, tdim)

30 boundary_facets = mesh.exterior_facet_indices(domain.topology)

31 boundary_dofs = fem.locate_dofs_topological(V, fdim,

32 boundary_facets)

33 bc = fem.dirichletbc(uD, boundary_dofs)

35 u = ufl.TrialFunction(V)

36 v = ufl.TestFunction(V)

38 ## termo de fonte

39 x = ufl.SpatialCoordinate(domain)

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

41 f = 100.0 * ufl.exp(-100 * ((x[0] - 0.5)**2 + (x[1] - 0.5)**2))

43 ## forma bilinear

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

46 ## forma linear

47 L = f * v * ufl.dx

49 ## problema linear

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

51 petsc_options_prefix="linear_problem_primal",

52 petsc_options={

53 "ksp_type": "cg",

54 "pc_type": "hypre",

55 "pc_hypre_type": "boomeramg",

56 "ksp_rtol": 1.0e-10,

57 },

58 )

59 uh = problem.solve()

61 # gráfico

62 import pyvista

63 from dolfinx import plot

65 # solução numérica

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

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

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

69 grid_uh.set_active_scalars("uh")

71 # plot

72 p = pyvista.Plotter()

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

74 p.add_axes()

75 p.add_legend()

76 p.view_xy()

77 p.show()

79 # dual problem

80 V_dual = fem.functionspace(domain, ('P', 3))

81 z = ufl.TrialFunction(V_dual)

82 w = ufl.TestFunction(V_dual)

84 ## condição de contorno

85 import numpy as np

86 uD = fem.Function(V_dual)

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

89 tdim = domain.topology.dim

90 fdim = tdim - 1

91 domain.topology.create_connectivity(fdim, tdim)

92 boundary_facets = mesh.exterior_facet_indices(domain.topology)

93 boundary_dofs = fem.locate_dofs_topological(V_dual, fdim,

94 boundary_facets)

95 bc = fem.dirichletbc(uD, boundary_dofs)

97 ## forma bilinear

98 a = ufl.dot(ufl.grad(w), ufl.grad(z)) * ufl.dx

100 ## funcional objetivo: média de u

101 M = w * ufl.dx

102

103 ## problema linear

104 problem = LinearProblem(a, M, bcs=[bc],

105 petsc_options_prefix="linear_problem_dual",

106 petsc_options={

107 "ksp_type": "cg",

108 "pc_type": "hypre",

109 "pc_hypre_type": "boomeramg",

110 "ksp_rtol": 1.0e-10,

111 })

112 zh = problem.solve()

113

114 # DWR indicator: eta_K ~ R_K(zh - Ih zh)

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

116 w0 = ufl.TestFunction(W0)

117

118 # Project dual solution to primal space to build the DWR weight

119 z_proj = fem.Function(V)

120 z_proj.interpolate(zh)

121 z_tilde = zh - z_proj

122

123 # Strong residual of -Delta(u) = f, weighted by dual error

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

125

126 # Flux jump term, also weighted by dual error

127 n = ufl.FacetNormal(domain)

128 flux_jump = ufl.jump(ufl.dot(ufl.grad(uh), n))

129 eta_form = fem.form(

130 r * z_tilde * w0 * ufl.dx

131 + 0.5 * flux_jump * ufl.avg(z_tilde) * (w0('+') + w0('-')) * ufl.dS

132 )

133

134 eta_vec = create_vector(fem.extract_function_spaces(eta_form))

135 with eta_vec.localForm() as loc_eta:

136 loc_eta.set(0)

137 assemble_vector(eta_vec, eta_form)

138 eta_vec.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,

139 mode=PETSc.ScatterMode.REVERSE)

140 eta_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES,

141 mode=PETSc.ScatterMode.FORWARD)

142

143 print(f"Refinement {i}: max |cell DWR indicator| = {np.abs(eta_vec.array).max()}")

144

145 # plot cell indicators (DG0 is cellwise, so use mesh topology + cell_data)

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

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

148 grid_residual.cell_data["cell_eta"] = np.abs(eta_vec.array.real)

149 grid_residual.set_active_scalars("cell_eta")

150

151 pl = pyvista.Plotter()

152 pl.add_mesh(grid_residual, show_edges=True, scalars="cell_eta", cmap="inferno", label="cell_eta")

153 pl.add_axes()

154 pl.add_legend()

155 pl.view_xy()

156 pl.show()

157

158 # mark cells for refinement based on residuals

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

160 tdim = domain.topology.dim

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

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

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

164 for c in local_cells:

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

166 cell_eta[c] = abs(eta_vec.array[dof].real)

167

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

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

170 threshold = 0.5 * global_max

171 cells_to_refine = local_cells[cell_eta > threshold]

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

173

174 # # global mesh refinement

175 # from dolfinx.mesh import refine

176 # refine_out = refine(domain)

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

178

179 # local mesh refinement

180 from dolfinx.mesh import refine

181 if len(cells_to_refine) > 0:

182 marker_entities = cells_to_refine

183 try:

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

185 marker_entities = mesh.compute_incident_entities(domain.topology,

186 cells_to_refine,

187 tdim, 1)

188 except Exception:

189 pass

190 refine_out = refine(domain, marker_entities)

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

192

193p.screenshot("fig.png")

194

195# crop to content and get size

196import os

197os.system("convert fig.png -trim fig.png")

198from PIL import Image

199im = Image.open("fig.png")

200w, h = im.size

201dpi = im.info.get('dpi', (300, 300))

202print(f'w={w/dpi[0]:.3f}in, h={h/dpi[1]:.3f}in')

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.

Política de uso de dados

Política de uso de dados

Método dos Elementos Finitos

2.4 Malha adaptativa dual

Teorema 2.4.1.(Estimativa a posteriori dual)

Envie seu comentário