Cubic mesh, problems with Refine and Assemble routines.

Sepulveda
Topic Author
Offline
New Member

7 years 2 months ago #156 by Sepulveda

Cubic mesh, problems with Refine and Assemble routines. was created by Sepulveda

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:

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.

joachim
Offline
Administrator

7 years 2 months ago #157 by joachim

Replied by joachim on topic Cubic mesh, problems with Refine and Assemble routines.

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

Code:

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))

ddrake
Offline
Premium Member

7 years 2 months ago #158 by ddrake

Replied by ddrake on topic Cubic mesh, problems with Refine and Assemble routines.

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.

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)

Best,
Dow

cwinters
Offline
Administrator

7 years 2 months ago #171 by cwinters

Replied by cwinters on topic Cubic mesh, problems with Refine and Assemble routines.

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".

Code:

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)

The following user(s) said Thank You: ddrake

Time to create page: 0.098 seconds