Cubic mesh, problems with Refine and Assemble routines.

Main Forum Old-Forum - do not post here
1

Hello,

I installed NGSolve last Monday with -DUSE_UMFPACK=ON, and I have these problems when running a 3D problem with a cubic mesh.

  1. I cannot use mesh.Refine()
  2. I have problems with a.Assemble, when using direct preconditioner

For example, running this code:

from ngsolve import *
from netgen.csg import unit_cube
mesh=Mesh(unit_cube.GenerateMesh(maxh=0.125, quad_dominated=True))
#mesh.Refine() #This is not working

H = H1(mesh, order=1)


u = H.TrialFunction()
v = H.TestFunction()

a = BilinearForm(H)
a+= SymbolicBFI( u*v)


c = Preconditioner(a, type="direct")

SetHeapSize(int(1e8))
a.Assemble()

It returns:

a.Assemble()
RuntimeError: caught exception in DirectPreconditioner:
UmfpackInverse: Numeric factorization failed.
needs a sparse matrix (or has memory problems)in Assemble BilinearForm ‘bfa’

Is any other way to make this work?
I also tried a.Assemble(heapsize=int(1e9)).

Thanks,
Paulina.

2

Hi Paulina,

Netgen cannot generate hex meshes. The mesher gives an error output, but the py-script continues until the solver complains about inverting the 0-matrix.

In the upcoming release we also catch Netgen exceptions such that one sees the origin of the problem directly.

If you want to generate a structured hex-mesh of a cube you can generate it by the py-script below,

best, Joachim

from netgen.meshing import *

def GenerateCubeMesh (n):
    
    mesh = Mesh()

    n = 10

    pids = []
    for i in range(n+1):
        for j in range(n+1):
            for k in range(n+1):
                pids.append (mesh.Add (MeshPoint(Pnt(i/n, j/n, k/n))))
            
    for i in range(n):
        for j in range(n):
            for k in range(n):
                base = i * (n+1)*(n+1) + j*(n+1) + k
                baseup = base+(n+1)*(n+1)
                pnum = [base,base+1,base+(n+1)+1,base+(n+1),
                        baseup, baseup+1, baseup+(n+1)+1, baseup+(n+1)]
                elpids = [pids[p] for p in pnum]
                mesh.Add (Element3D(1,elpids))

    def AddSurfEls (p1, dx, dy, facenr):
        for i in range(n):
            for j in range(n):
                base = p1 + i*dx+j*dy
                pnum = [base, base+dx, base+dx+dy, base+dy]
                elpids = [pids[p] for p in pnum]
                mesh.Add (Element2D(facenr,elpids))
                        
    mesh.Add (FaceDescriptor(surfnr=1,domin=1,bc=1))
    AddSurfEls (0, 1, n+1, facenr=1)

    mesh.Add (FaceDescriptor(surfnr=2,domin=1,bc=2))
    AddSurfEls (0, (n+1)*(n+1), 1, facenr=2)

    mesh.Add (FaceDescriptor(surfnr=3,domin=1,bc=3))
    AddSurfEls (0, n+1, (n+1)*(n+1), facenr=3)

    mesh.Add (FaceDescriptor(surfnr=4,domin=1,bc=4))
    AddSurfEls ((n+1)**3-1, -(n+1), -1, facenr=1)

    mesh.Add (FaceDescriptor(surfnr=5,domin=1,bc=5))
    AddSurfEls ((n+1)**3-1, -(n+1)*(n+1), -(n+1), facenr=1)

    mesh.Add (FaceDescriptor(surfnr=6,domin=1,bc=6))
    AddSurfEls ((n+1)**3-1, -1, -(n+1)*(n+1), facenr=1)

    return mesh




mesh = GenerateCubeMesh(10)
mesh.Save ("cubemesh.vol")

import ngsolve
ngsolve.Draw (ngsolve.Mesh(mesh))
3

The code works to generate a quad-dominated unit cube, but could there be a typo in setting the indices for the last three arrays of faces?

The first three calls to AddSurfEls() are passing sequential values for facenr, matching the surfnr of the preceding FaceDescriptor() call, but the last three are passing the value 1.

[quote][code]
mesh.Add (FaceDescriptor(surfnr=1,domin=1,bc=1))
AddSurfEls (0, 1, n+1, facenr=1)

mesh.Add (FaceDescriptor(surfnr=2,domin=1,bc=2))
AddSurfEls (0, (n+1)*(n+1), 1, facenr=2)

mesh.Add (FaceDescriptor(surfnr=3,domin=1,bc=3))
AddSurfEls (0, n+1, (n+1)*(n+1), facenr=3)

mesh.Add (FaceDescriptor(surfnr=4,domin=1,bc=4))
AddSurfEls ((n+1)**3-1, -(n+1), -1, facenr=1)

mesh.Add (FaceDescriptor(surfnr=5,domin=1,bc=5))
AddSurfEls ((n+1)**3-1, -(n+1)*(n+1), -(n+1), facenr=1)

mesh.Add (FaceDescriptor(surfnr=6,domin=1,bc=6))
AddSurfEls ((n+1)**3-1, -1, -(n+1)*(n+1), facenr=1)

[/code][/quote]
Best,
Dow

4

You are right, there shouldn’t be a facenr=1 at the last three calls of AddSurfEls(). This lead to a mesh with three boundaries with the numbers 1,2 and 3. The last three FaceDescriptors were unused.

The “facenr” has to match the number of the FaceDescriptor. If you call mesh.Add() with a FaceDescriptor as argument, you get it’s number as return value. Meaning you could use the return value as “facenr”.

fd1 = mesh.Add (FaceDescriptor(surfnr=1,domin=1,bc=1))
AddSurfEls (0, 1, n+1, facenr=fd1)

fd2 = mesh.Add (FaceDescriptor(surfnr=2,domin=1,bc=2))
AddSurfEls (0, (n+1)*(n+1), 1, facenr=fd2)

fd3 = mesh.Add (FaceDescriptor(surfnr=3,domin=1,bc=3))
AddSurfEls (0, n+1, (n+1)*(n+1), facenr=fd3)

fd4 = mesh.Add (FaceDescriptor(surfnr=4,domin=1,bc=4))
AddSurfEls ((n+1)**3-1, -(n+1), -1, facenr=fd4)

fd5 = mesh.Add (FaceDescriptor(surfnr=5,domin=1,bc=5))
AddSurfEls ((n+1)**3-1, -(n+1)*(n+1), -(n+1), facenr=fd5)

fd6 = mesh.Add (FaceDescriptor(surfnr=6,domin=1,bc=6))
AddSurfEls ((n+1)**3-1, -1, -(n+1)*(n+1), facenr=fd6)