Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Cubic mesh, problems with Refine and Assemble routines.

More
6 years 7 months ago #156 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.
More
6 years 7 months ago #157 by joachim
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))
More
6 years 7 months ago #158 by ddrake
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
More
6 years 6 months ago #171 by cwinters
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.166 seconds