firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
512 stars 160 forks source link

Labelling of Internal Facets in Mesh Generated with Netgen #3668

Open adofarsi opened 4 months ago

adofarsi commented 4 months ago

Describe the current issue Currently, when generating meshes using Netgen, there is no functionality available to correctly label the internal facets. This lack of labelling can lead to difficulties in identifying and differentiating between various internal facets of the mesh for DG discretisations.

Describe the solution you'd like Addition of a feature that enables the correct labelling of internal facets in the mesh generated by Netgen.

Additional info Code for testing below:

import netgen.gui
from netgen.occ import * 
from firedrake import *

outer = Rectangle(1, 1).Face()
outer.edges.name="outer"
outer.edges.Max(X).name = "r"
outer.edges.Min(X).name = "l"
outer.edges.Min(Y).name = "b"
outer.edges.Max(Y).name = "t"

Inner = MoveTo(0.1, 0.1).Rectangle(0.3, 0.5).Face()
Inner.edges.name="interface"
outer = outer - Inner

Inner.faces.name="inner"
Inner.faces.col = (1, 0, 0)
outer.faces.name="outer"

geo = Glue([Inner, outer])
ngmesh = OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.2)
mesh = Mesh(ngmesh)

V = FunctionSpace(mesh, "DG", 2)
u = Function(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = Function(V).interpolate(sin(x[0]*pi)*sin(2*x[1]*pi))
pen = Constant(1e10)
F = + inner(grad(u), grad(v)) * dx \
    - f * v * dx \
    + pen * (u-2) * v * ds(8) \
    + pen * jump(u) * jump(v) * dS(1) \
    + pen * jump(u) * jump(v) * dS(2) \
    + pen * jump(u) * jump(v) * dS(3) \
    + pen * jump(u) * jump(v) * dS(4) \
    + pen * jump(u) * jump(v) * dS(5) \
    + pen * jump(u) * jump(v) * dS(6) 
bcs = [DirichletBC(V, Constant(2.0), (1,))]
solve(F == 0, u, bcs=[], solver_parameters={"ksp_type": "preonly",
                                            "pc_type": "lu"})
File("testing.pvd").write(u)
connorjward commented 3 months ago

@UZerbinati how easy do you think this will be to do?

UZerbinati commented 3 months ago

@connorjward Is a work in progress in this branch: https://github.com/NGSolve/ngsPETSc/tree/uz/dglabels. It works for 2D but now needs to be extended to 3D. The main issue is that it changes the labelling routine we have used so far. I think it could be useful to create a labelling utility function directly in ngsPETSc that retrieves from a boundary string the corresponding ID. Maybe if a mesh is a netgen mesh we can use this utility function directly in Firdreke DirichletBC?

connorjward commented 3 months ago

@connorjward Is a work in progress in this branch: https://github.com/NGSolve/ngsPETSc/tree/uz/dglabels. It works for 2D but now needs to be extended to 3D. The main issue is that it changes the labelling routine we have used so far. I think it could be useful to create a labelling utility function directly in ngsPETSc that retrieves from a boundary string the corresponding ID. Maybe if a mesh is a netgen mesh we can use this utility function directly in Firdreke DirichletBC?

Maybe. Does DMLabel not support this kind of thing? I think we want to keep netgen code localised to mesh.py where possible.