How to compose a region using a gfu (or cf)?

Main Forum Old-Forum - do not post here
1

Hello,

I don’t understand how i could compose a ngsolve.comp.Region using a given gfu (or cf). Using that region i’d like to separate two groups of dofs (using the definedon=region).

Here the concrete goal:

  • sort every dof of fes into a second fes where the gfu value is higher than tol. All the other dofs should go into a third fes.

Perhaps its easy and I just didn’t go the correct way of solving the problem.
Thanks for your help!

2

Hi creativeworker,

you can use a BitArray to sort dofs. This works only for lowest order H1 elements where every dof is associated with the vertex value.

For example

fes = H1(mesh, order=1)
gfu = GridFunction(fes)
gfu.Set(x*y+2*y)

limit = 1.5

dofs = BitArray(fes.ndof)
dofs[:] = False
for i in range(fes.ndof):
    if gfu.vec[i] > limit:
        dofs[i] = True

saves all dofs where gfu is greater than 1.5.

Defining a second space which lives only on these dofs can be done only with

fes2 = H1(mesh, order=1, dirichlet="left|top|bottom|right")
for i in range(fes.ndof):
    if dofs[i] == False:
        fes2.SetCouplingType(i,COUPLING_TYPE.UNUSED_DOF)
fes2 = Compress(fes2)

Best,
Michael

https://ngsolve.org/media/kunena/attachments/889/sortdof.py

Attachment: sortdof_2020-04-22-2.py

3

Hi,

Short addition. Compress has an argument “active_dofs” exactly for replacing the last loop:

fes2 = H1(mesh, order=1, dirichlet="left|top|bottom|right")
fes2 = Compress(fes2, active_dofs=dofs)

Best,
Christoph

4

Thank you very much. The Compress function helps me also in another place!

5

May i ask two further things on that topic:

Is there a possibility to find the boundary dofs of such a bit mask? Or to define a region with the bit mask and find the boundary dofs there?

Secondly, is there a chance to bring this idea to HCurl space?

6

Hi creativeworker,

you can get the dofs of a boundary as a BitArray with

fes.GetDofs(mesh.Boundaries("left"))

which then can be combined with other BitArrays, see attached file.

Best
Michael

https://ngsolve.org/media/kunena/attachments/889/sortdof_bnd.py

Attachment: sortdof_bnd.py

7

Thanks for the fast reply.
What I meant was more the other way round. So I have created my bit mask (using the limit). Now I need the dofs that surround the areas of where the limit of the gfu was beyond the limit speaking in the physical space. So, the dofs have to somehow be interpreted back in the mesh domain and decided wether the dof has only marked dofs around itself or not.

8

So, if I understood you correctly you would like to have some kind of level-set function.

You can write a code finding the vertices on the boundary, this tutorial might help looping over elements, edges, etc.

Or you might have a look into the ngsxfem package, where a level-set can be defined (but there elements may be cut instead of using the vertices, but I’m not an expert for this package)

9

Thank you! I tried with iterating over all the elements. But then something strange happens. The second for-loop is stuck on the last element. So the elements are not iterated from start to end again.

Here is my code-snippet:

        elements = list(fes.Elements())
        masterelements = []

        for i in range(fes.ndof):
            if gfu.vec[i] > limit:
                for el in elements:
                    if i in el.dofs:
                        masterelements.append(el)
10

can you post a complete failing example ?
Joachim

11

The attached example shows my problem.
https://ngsolve.org/media/kunena/attachments/934/element-iteration_min.py

Attachment: element-iteration_min.py

12

Making a python-list from “fes.Elements()”

elements = list(fes.Elements())

does not work as you would expect it to. The reason is that fes.Elements() is only meant to be iterated through once, and is not a real “array” on C++ side.

This is not a good way to organize the loop - you go through every element for every DOF. Also, some elements will be added multiple times. Also, again, making lists from elements does not work, they are only meant to be temporary, not persistent.

for i in range(fes.ndof):
    if gfu.vec[i] > limit:
        for el in elements:
            if i in el.dofs:
                masterelements.append(el)

Better to go through elements and for each element to only check it’s DOFs, and only save the numbers of elements.

master_el_nrs = []
for el in fes.Elements():
    is_master = False
    for dof in el.dofs:
        is_master |= gfu.vec[dof] > limit
    if is_master:
        master_el_nrs.append(el.nr)

We can always iterate through the numbers, and then get the elements with those numbers:

# Iterate through chosen elements:
for el_nr in master_el_nrs:
    el_id = ElementId(VOL, el_nr)
    print("VOL Element #", el_nr, "was chosen, dofs =", fes.GetDofNrs(el_id))
    el = fes.GetFE(el_id)
    # now we could do something with the element...

Hope this helps.

Best,
Lukas

Attachment: element-iteration_min.py

13

With the coming nightly version you can keep the FESpace-elements in a list, and obtain dof numbers from them.

I am completely in favour of Lukas’ O(N) algorithm (instead of the original O(N^2) method). With the nightly fixes, you can store numbers, element-ids, or FESpace-elements in the list.

Joachim