add delaunay mesh region on the top of prism(slice) surface

asilalahi · January 11, 2019, 3:50am

Hi,

I was able to generate a 3-region cylindrical shells with parameters:

490<r<500 , Region 1
500<r<505 , Region 2 (membrane)
505<r<515 , Region 3

Here, the Region 1 and 3 are separated by a 5nm width membrane.

The code is as bellow :

from netgen.csg import *

geo = CSGeometry()
cube = OrthoBrick(Pnt(-5700,-5700,0), Pnt(7100,7100,6300))

cyl1_in=Cylinder(Pnt(700,700,0), Pnt(700,700,200),490)
cyl1=Cylinder(Pnt(700,700,0), Pnt(700,700,200),500)
cyl2 = Cylinder(Pnt(700,700,0), Pnt(700,700,200),505)
cyl3 = Cylinder(Pnt(700,700,0), Pnt(700,700,200),515)

dom1=cyl1 * cube - cyl1_in
dom2=cyl2 * cube - cyl1
dom3=cyl3 * cube - cyl2

geo.Add(dom1)
geo.Add(dom2)
geo.Add(dom3)

n =10
slices1x = [i/n for i in range(1,n)]

slices2 = [i/2 for i in range(1,2)]

geo.CloseSurfaces(cyl1,cyl1_in,slices1x)
geo.CloseSurfaces(cyl1,cyl2,slices2)
geo.CloseSurfaces(cyl2,cyl3,slices1x)

geo.SetBoundingBox(Pnt (-12000, -12000, -12000), Pnt(12000, 12000, 12000))

ngmesh = geo.GenerateMesh(maxh=100)
ZRefinement(ngmesh, geo)

from ngsolve import *
mesh = Mesh(ngmesh)
Draw(mesh)

===============================================================================

The most important & signifficant physics occur in these 3 regions, but i still need to put correct boundary, especially outside Region 3.
For this purpose, I would like to add another region/domain using delaunay mesh outside Region 3 , maybe up to r=3000nm. How to do this without my mesh being constrained by maxh=100( as in above code) ? I would like to make it coarse in this region, maybe maxh=700.

Best,
Alex

cwinters · January 11, 2019, 8:27am

Hi,

you can define the maxh right when you add your domains to the geometry

geo.Add(dom1,maxh=100)
geo.Add(dom2,maxh=100)
geo.Add(dom3,maxh=100)

If you add another cylinder on the outside without adding CloseSurfaces you get a tetrahedral mesh in this domain.

ngmesh = geo.GenerateMesh(maxh=700)

should then give you a courser mesh on the outside domain.

Best,
Christoph

Attachment: regions.py

asilalahi · January 11, 2019, 9:18am

Hi Christopher,

Thanks alot, i didnt know that.

I have another setup that i want to do, in which i want to add a pillar (cylindrical geometry) to the previous geometry. This pillar is located right outside Region 2. And because we also have region 3 outside Region 2, there will be overlapping between Region 3 and pillar.

So here is the code that i tried:

===================================================================================
from netgen.csg import *

geo = CSGeometry()
cube = OrthoBrick(Pnt(-5700,-5700,0), Pnt(7100,7100,1300))

cube_b = OrthoBrick(Pnt(-5700,-100,0), Pnt(7100,195,2000))

cyl1_in=Cylinder(Pnt(700,700,0), Pnt(700,700,200),490)
cyl1=Cylinder(Pnt(700,700,0), Pnt(700,700,200),500)
cyl2 = Cylinder(Pnt(700,700,0), Pnt(700,700,200),505)
cyl3 = Cylinder(Pnt(700,700,0), Pnt(700,700,200),550)

cyl4=Cylinder(Pnt(700,0,650), Pnt(700,-700,650),200)

dom1=cyl1 * cube - cyl1_in
dom2=cyl2 * cube - cyl1
dom3=cyl3 * cube - cyl2

dom4=cyl4* cube_b

geo.Add(dom1.maxh(100))
geo.Add(dom2.maxh(100))
geo.Add(dom3.maxh(100))
geo.Add(dom4.maxh(50))

n =10
slices1x = [i/n for i in range(1,n)]

slices2 = [i/2 for i in range(1,2)]

geo.CloseSurfaces(cyl1,cyl1_in,slices1x)
geo.CloseSurfaces(cyl1,cyl2,slices2)
geo.CloseSurfaces(cyl2,cyl3,slices1x)

geo.SetBoundingBox(Pnt (-12000, -12000, -12000), Pnt(12000, 12000, 12000))

ngmesh = geo.GenerateMesh(maxh=100)
ZRefinement(ngmesh, geo)

from ngsolve import *
mesh = Mesh(ngmesh)
Draw(mesh)

==============================================================================
the geometry roughly looks like :

Obviously Region 3 (exterior) and Region 4( the pillar) overlap, as shown below:

I want to first try the easier way to get around this problem, which is to delete region 3 that intersects with region 4, while maintaining the parallel prism cell in region 3. How to do this?

Thanks,
Alex

cwinters · January 11, 2019, 12:24pm

Hi,

you have two tell the mesher how the two parts are “connected”. In your case you have to cut out “cyl3”:

dom4=cyl4* cube_b - cyl3

Best,
Christoph

cwinters · January 11, 2019, 12:46pm

Hi,

I forgot to add that you can not connect your pillar (cyl4) to cyl3 when the volume between cyl2 and cyl3 consists of prisms. The intersection of the cylinders has to be resolved in the surface mesh.To add the prisms between two layers, the meshes on these layer must be identical. This is not the case if one intersects with another cylinder while the other one does not.

Best,
Christoph

asilalahi · January 11, 2019, 1:12pm

Hi Christopher,

Thanks for your reply, yes it is the challenge with my setup.
To resolve the important physics around membrane boundary, i need to have resolution around ~ 0.5 nm in radial direction near the membrane. Thus if i dont use prism cell, i will need too many number of tetrahedral cells in Region 3. because region 3 range from 505 to 515 nm and with length of 1300nm.

I thought that if i could delete some of the prism cells in region 3 that intersect with pillar, then it maybe good enough for my first try.

FYI:
Region 2: membrane (physical)
Region 4: pillar (pyhsical)

Region 1 and 3 represent mixture of ions that can freely move.

Best,
Alex