Hi everyone,

I am trying to get a local solution from a DG implementation. To do that, I require only the inner information (’+’) in the interior facet integration (i.e. dS) instead of using the jumps and averages.

The issue is, apparently, the facet normal directions (’+’) change when I use the function refine(), but not when I do a manual uniform refinement (creating a new refined mesh on each level).

To illustrate, I wrote this short code, where I calculate the L2 norm using 1) a manual uniform refinement and 2) using the refine() function.

For the first case I have a convergent behaviour in the L2 norm

l2norm: [19.0606515, 1.6594006, 0.5111520, 0.37043388, 0.331217411, 0.316423740].

But in the second case (same code just commenting the line 37), I have the results I want only for the first two solutions (first refinement), because then the global matrix is not the same anymore.

l2norm: [19.06065, 1.6594006, inf, nan, inf, inf].

If you have any suggestions on how to fix the problem, I would appreciate it!

Thank you,

Juan

CODE:

```
1. import fenics
2. from dolfin import *
3. import numpy as np
4. def bh1(u,v):
5. val = inner(kappa*grad(u), grad(v))*dx
6. val += -inner(u, dot(kappa*grad(v), n))*ds
7. val += -inner(v, dot(kappa*grad(u), n))*ds
8. val += inner((eta/h*kappa)*v, u)*ds
9. val += dot(beta_vec, grad(u))*v*dx
10. val += 0.5*(abs(dot(beta_vec,n))-dot(beta_vec,n))*u*v*ds
11. return val
12. def bhfacet(u,v):
13. val = -inner(dot(avg(kappa*grad(u)), n('+')), jump(v))*dS
14. val += -inner(jump(u), dot(avg(kappa*grad(v)), n('+')))*dS
15. val += inner((avg(eta/h)*avg(kappa))*jump(u), jump(v))*dS
16. val -= dot(beta_vec('+'),n('+'))*jump(u)*avg(v)*dS
17. val += 0.5*abs(dot(beta_vec('+'),n('+')))*jump(u)*jump(v)*dS
18. return val
19. def bhfacetmod(u,v):
20. val = -inner(dot(kappa('+')*grad(u)('+'), n('+')), v('+'))*dS
21. val += -inner(u('+'), dot(kappa('+')*grad(v)('+'), n('+')))*dS
22. val += inner((eta/h('+')*kappa('+'))*u('+'), v('+'))*dS
23. val -= dot(beta_vec('+'),n('+'))*u('+')*v('+')*dS
24. val += 0.5*abs(dot(beta_vec('+'),n('+')))*u('+')*v('+')*dS
25. return val
26. def bhdg(u,v):
27. return bh1(u,v) + bhfacet(u,v)
28. def bhdglocal(u,v):
29. return bh1(u,v) + bhfacetmod(u,v)
30. kappa = Constant(1e-2) ; beta_vec= as_vector([1,]); eta = 18
31. nlist = [40,80,160,320,640,1280]
32. l2norm = []
33. mesh = IntervalMesh(nlist[0], 0, 1)
34. u_exa = Expression('(exp(1/1e-2*x[0])-1)/(exp(1/1e-2)-1)-x[0]',degree=5)
35. MAX_ITER = np.size(nlist)
36. for level in range(MAX_ITER):
37. mesh = IntervalMesh(nlist[level], 0, 1)
38. V = FunctionSpace(mesh, 'DG', 1)
39.
40. u = TestFunction(V)
41. v = TrialFunction(V)
42.
43. boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
44. boundary_markers.set_all(0)
45. ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
46.
47. n = FacetNormal(mesh); h = CellDiameter(mesh)
48. L = inner(v, -Constant(1))*dx
49.
50. solDG = Function(V)
51. solve(bhdg(u,v) == L ,solDG )
52.
53. bdrhs= bhdglocal(solDG,v)
54.
55. uDGmod = Function(V)
56. solve(bhdglocal(u,v) == L + bdrhs ,uDGmod)
57. l2norm.append(assemble(inner(u_exa-uDGmod,u_exa-uDGmod)*dx))
58. mesh = refine(mesh)
59.
60.
61. print('l2norm:',l2norm)
```