question on 3D adaptive mesh refinement

[Guosheng Fu]

Hello,

I have a question on 3D mesh refinement.
I manually mark a few elements for refinement, but the code seems just give me a uniform refinement. The following code works in 2D, giving me an adaptive mesh, but not in 3D.
I don't understand it....

from netgen.geom2d import unit_square
from netgen.csg import unit_cube
from time import sleep

ngsglobals.msg_level = 0

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.4))
#mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

l = []

def Mark():
    for el in mesh.Elements():    
      mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements():
        if el.nr < 5:
          mesh.SetRefinementFlag(el, True)

Draw(mesh)
input("yes?")

while mesh.ne < 100000:  
    Mark()
    mesh.Refine()
    Redraw()
    print(mesh.nv, mesh.ne)
    sleep(2)

Best,
Guosheng

[Joachim Schoeberl]

you can set refinement-flags for volume and boundary elements, you have to clear them both:

---

def Mark():
for el in mesh.Elements():
mesh.SetRefinementFlag(el, False)
for el in mesh.Elements(BND):
mesh.SetRefinementFlag(el, False)
for el in mesh.Elements():
if el.nr < 5:
mesh.SetRefinementFlag(el, True)

---

[Guosheng Fu]

Hi Joachim,

That works now. Thanks!

Here is another issue:
I want to use adaptive mesh refinement for HDG in 3D, but the FacetFEspace gives me a segmentation fault when updated on an adaptive mesh, which is fine on uniform refined mesh....
Here is the code that generates the segmentation fault.
I am quite interested in debugging the code in python, but I am a novince in debugging a code... Do you know how to backtrack the source of the bug?

from netgen.csg import unit_cube
from ngsolve import *
from time import sleep

ngsglobals.msg_level = 0

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.4))
fes = FacetFESpace(mesh=mesh, order = 1)

l = []

def Mark():
    for el in mesh.Elements():    
      mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements(BND):    
      mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements():
        if el.nr < 5:
          mesh.SetRefinementFlag(el, True)

Draw(mesh)

while mesh.ne < 100000:  
    Mark()
    mesh.Refine()
    fes.Update()
    Redraw()
    print(mesh.nv, mesh.ne)
    sleep(1)

Best,
Guosheng

[Joachim Schoeberl]

Hi Guosheng,

it was a bug, it's fixed with this commit:
https://gitlab.asc.tuwien.ac.at/jschoeberl/ngsolve/commit/b21b20561ca715b622bfcc600078cb5f6492d298

Joachim

[Guosheng Fu]

Hi Joachim,

Thanks! It works now.

Just another side question. I noticed there are quite a few "testout" messages in the source code helpful for debugging.
How can I print out those test out message in python?
I guess I shall compile ngsolve in debug mode, right? But I don't know how to do so..

Best,
Guosheng

[Joachim Schoeberl]

Hi Guosheng,

use in the py-script:

SetTestoutFile (filename)

secret note: You have to do it after mesh-generation, otherwise it's overwritten.

You can control the output via flags:
Add "-print" to the FESpace, or "-printelmat" to the BilinearForm flags.

Joachim

[Guosheng Fu]

Aha, I will never figure out these tricks myself ;>
Thanks!

[Sander Rhebergen]

Hi,
I have a related question. I would like to create on the unit square a 2D grid but refined near the boundary (to capture boundary layer effects) and coarser in the interior. For example, if I want a refined grid for all x < 0.1 then I thought I could use a slightly modified version of what Guosheng wrote above:

def Mark():
    for el in mesh.Elements(): 
        mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements(BND): 
        mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements():
        for v in el.vertices:
            if (v.x[0] < 0.1):
                mesh.SetRefinementFlag(el, True)

I understand that v.x[0] doesn't exist, but is there something available within ngsolve so that I can do something like this?
Is this the best way to create a grid with a boundary layer within NGSolve-Python?

[Joachim Schoeberl]

Hi Sander,

one way is to flag vertices on the boundary first, and then refine elements having the boundary-flag set. You get the boundary vertices by looping over boundary (n-1 dimensional elements), and loop over their vertices. Or you just use the fes.FreeDofs() bitarray of non-Dirichlet dofs for that.

But, some day we want to access point coordinates. NGSolve is missing that at the moment, but Netgen has it. You may have realized that there is a ngsove.Mesh, and a netgen.Mesh as well. GenerateMesh generates a netgen.Mesh, which we usually wrap into a ngsolve.Mesh (which has a shared pointer of the original netgen.Mesh). We can extract the netgen.Mesh, and get the point coordinates from it. Unfortunately, for historical reasons, the netgen-points ave 1-based indexing. You can do it like this.

from ngsolve import *
from netgen.geom2d import unit_square
from netgen.meshing import PointId

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

for i in range(mesh.nv):
    pnt = mesh.ngmesh[PointId(i+1)]
    print (pnt[0], pnt[1])

There is also the possibility to use anisotropic refinement for boundary layers in Netgen/NGSolve, more follows,

Joachim

[Sander Rhebergen]

Hi Joachim,
Thanks, that was what I needed. Hereby the code that creates a mesh with refinement near the boundary:

from ngsolve import *
from netgen.geom2d import unit_square
from netgen.meshing import PointId

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

bl = 0.02

def Mark():
    for el in mesh.Elements(): 
        mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements(BND): 
        mesh.SetRefinementFlag(el, False)
    for el in mesh.Elements():
        for i in el.vertices:
            pnt = mesh.ngmesh[PointId(i+1)]
            if (pnt[0] < bl or pnt[0] > 1-bl or pnt[1] < bl or pnt[1] > 1-bl):
                mesh.SetRefinementFlag(el, True)

Draw(mesh)    

while mesh.ne < 10000:  
    Mark()
    mesh.Refine()
    Redraw()
    print(mesh.nv, mesh.ne)

Sander