from ngsolve import *
from netgen.csg import *
from ngsolve.internal import *


def ATCerror_brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    box = OrthoBrick(Pnt(-Rminus,-Rminus,-Rminus),Pnt(Rminus,Rminus,Rminus+delta))
    pl1 = Plane(Pnt(0,0,Rminus),Vec(0,0,-1)).bc("crack")
    pl2 = Plane(Pnt(0,0,Rminus+delta),Vec(0,0,1))#.bc("top")
    o_minus = (box - pl1)
    geometry.Add ((box - pl1).mat("ominus"),bcmod=[(o_minus,"nocrack")])    
    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-box).mat("oplus").maxh(hmax))
    geometry.Add ((box * pl1 * pl2).mat("olayer").maxh(hmax),bcmod=[(pl1,"crack"),(box,"sides"),(pl2,"top")])
    #slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2)#,slices)
    
    return geometry

def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus+.2).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    circle = Sphere(Pnt(0,0,0),Rminus)
    small_cylinder = Cylinder(Pnt(0,0,1),Pnt(0,0,-1),Rminus*0.5)

    pl1 = Plane(Pnt(0,0,0),Vec(0,0,-1))
    pl2 = Plane(Pnt(0,0,delta),Vec(0,0,1))#.bc("top")

    o_minus_with_layer = (circle*pl2)
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = ((pl1 * pl2 * circle) - small_cylinder)
    yes_olayer =(pl1 * pl2 * small_cylinder)

    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-o_minus_with_layer).mat("oplus").maxh(hmax))
    geometry.Add ((o_minus).mat("ominus").maxh(hmax))    
    geometry.Add (yes_olayer.mat("olayer").maxh(hmax),bcmod=[(pl1,"crack")])
    geometry.Add (no_olayer.mat("oplus").maxh(hmax))
    slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2,slices)
    
    return geometry

def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")

    pml = Sphere(Pnt(0,0,0),Rpml)

    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    box = OrthoBrick(Pnt(-Rminus,-Rminus,-Rminus),Pnt(Rminus,Rminus,Rminus))
    pl1 = Plane(Pnt(0,0,Rminus),Vec(0,0,-1)).bc("crack")
    pl2 = Plane(Pnt(0,0,Rminus+delta),Vec(0,0,1))#.bc("top")
    o_minus = (box - pl1).maxh(hsample)

    geometry.Add (o_minus.mat("ominus").maxh(hmax))    
    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-box).mat("oplus").maxh(hmax))
    geometry.Add ((box * pl1 * pl2).mat("olayer").maxh(hmax))

    #slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2)#,slices)
    
    return geometry

def sphere_geometry(Rminus, Rplus, Rext, Rpml, Rother, c, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    
    o_minus = (Sphere(Pnt(0,0,0), Rminus)).maxh(hmax).mat("ominus").maxh(hsample)
    other = Sphere(Pnt(0,c,0),Rother)
    withCrack = (o_minus * other)
    withoutCrack = (o_minus - other)

    geometry.Add ((o_ext - pml).mat("air"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-o_minus-other).mat("oplus"))
    geometry.Add ((other-o_minus).mat("olayer"))
    geometry.Add (withCrack,bcmod=[(o_minus,"crack")])
    geometry.Add (withoutCrack,bcmod=[(o_minus,"nocrack")])

    return geometry

b528d361a5611208e4b9da92b76342bb2d9b420c
<mask token>


def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus + 0.2).bc('interface')
    circle = Sphere(Pnt(0, 0, 0), Rminus)
    small_cylinder = Cylinder(Pnt(0, 0, 1), Pnt(0, 0, -1), Rminus * 0.5)
    pl1 = Plane(Pnt(0, 0, 0), Vec(0, 0, -1))
    pl2 = Plane(Pnt(0, 0, delta), Vec(0, 0, 1))
    o_minus_with_layer = circle * pl2
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = pl1 * pl2 * circle - small_cylinder
    yes_olayer = pl1 * pl2 * small_cylinder
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus_with_layer).mat('oplus').maxh(hmax))
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add(yes_olayer.mat('olayer').maxh(hmax), bcmod=[(pl1, 'crack')])
    geometry.Add(no_olayer.mat('oplus').maxh(hmax))
    slices = [(2 ** -i) for i in reversed(range(1, 6))]
    geometry.CloseSurfaces(pl1, pl2, slices)
    return geometry


def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus,
        Rminus))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = (box - pl1).maxh(hsample)
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax))
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


<mask token>

<mask token>


def ATCerror_brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus, 
        Rminus + delta))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = box - pl1
    geometry.Add((box - pl1).mat('ominus'), bcmod=[(o_minus, 'nocrack')])
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax), bcmod=[(pl1,
        'crack'), (box, 'sides'), (pl2, 'top')])
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus + 0.2).bc('interface')
    circle = Sphere(Pnt(0, 0, 0), Rminus)
    small_cylinder = Cylinder(Pnt(0, 0, 1), Pnt(0, 0, -1), Rminus * 0.5)
    pl1 = Plane(Pnt(0, 0, 0), Vec(0, 0, -1))
    pl2 = Plane(Pnt(0, 0, delta), Vec(0, 0, 1))
    o_minus_with_layer = circle * pl2
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = pl1 * pl2 * circle - small_cylinder
    yes_olayer = pl1 * pl2 * small_cylinder
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus_with_layer).mat('oplus').maxh(hmax))
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add(yes_olayer.mat('olayer').maxh(hmax), bcmod=[(pl1, 'crack')])
    geometry.Add(no_olayer.mat('oplus').maxh(hmax))
    slices = [(2 ** -i) for i in reversed(range(1, 6))]
    geometry.CloseSurfaces(pl1, pl2, slices)
    return geometry


def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus,
        Rminus))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = (box - pl1).maxh(hsample)
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax))
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


<mask token>

<mask token>


def ATCerror_brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus, 
        Rminus + delta))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = box - pl1
    geometry.Add((box - pl1).mat('ominus'), bcmod=[(o_minus, 'nocrack')])
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax), bcmod=[(pl1,
        'crack'), (box, 'sides'), (pl2, 'top')])
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus + 0.2).bc('interface')
    circle = Sphere(Pnt(0, 0, 0), Rminus)
    small_cylinder = Cylinder(Pnt(0, 0, 1), Pnt(0, 0, -1), Rminus * 0.5)
    pl1 = Plane(Pnt(0, 0, 0), Vec(0, 0, -1))
    pl2 = Plane(Pnt(0, 0, delta), Vec(0, 0, 1))
    o_minus_with_layer = circle * pl2
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = pl1 * pl2 * circle - small_cylinder
    yes_olayer = pl1 * pl2 * small_cylinder
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus_with_layer).mat('oplus').maxh(hmax))
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add(yes_olayer.mat('olayer').maxh(hmax), bcmod=[(pl1, 'crack')])
    geometry.Add(no_olayer.mat('oplus').maxh(hmax))
    slices = [(2 ** -i) for i in reversed(range(1, 6))]
    geometry.CloseSurfaces(pl1, pl2, slices)
    return geometry


def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus,
        Rminus))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = (box - pl1).maxh(hsample)
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax))
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


def sphere_geometry(Rminus, Rplus, Rext, Rpml, Rother, c, delta, hsample, hmax
    ):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    o_minus = Sphere(Pnt(0, 0, 0), Rminus).maxh(hmax).mat('ominus').maxh(
        hsample)
    other = Sphere(Pnt(0, c, 0), Rother)
    withCrack = o_minus * other
    withoutCrack = o_minus - other
    geometry.Add((o_ext - pml).mat('air'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus - other).mat('oplus'))
    geometry.Add((other - o_minus).mat('olayer'))
    geometry.Add(withCrack, bcmod=[(o_minus, 'crack')])
    geometry.Add(withoutCrack, bcmod=[(o_minus, 'nocrack')])
    return geometry

from ngsolve import *
from netgen.csg import *
from ngsolve.internal import *


def ATCerror_brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus, 
        Rminus + delta))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = box - pl1
    geometry.Add((box - pl1).mat('ominus'), bcmod=[(o_minus, 'nocrack')])
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax), bcmod=[(pl1,
        'crack'), (box, 'sides'), (pl2, 'top')])
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus + 0.2).bc('interface')
    circle = Sphere(Pnt(0, 0, 0), Rminus)
    small_cylinder = Cylinder(Pnt(0, 0, 1), Pnt(0, 0, -1), Rminus * 0.5)
    pl1 = Plane(Pnt(0, 0, 0), Vec(0, 0, -1))
    pl2 = Plane(Pnt(0, 0, delta), Vec(0, 0, 1))
    o_minus_with_layer = circle * pl2
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = pl1 * pl2 * circle - small_cylinder
    yes_olayer = pl1 * pl2 * small_cylinder
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus_with_layer).mat('oplus').maxh(hmax))
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add(yes_olayer.mat('olayer').maxh(hmax), bcmod=[(pl1, 'crack')])
    geometry.Add(no_olayer.mat('oplus').maxh(hmax))
    slices = [(2 ** -i) for i in reversed(range(1, 6))]
    geometry.CloseSurfaces(pl1, pl2, slices)
    return geometry


def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    box = OrthoBrick(Pnt(-Rminus, -Rminus, -Rminus), Pnt(Rminus, Rminus,
        Rminus))
    pl1 = Plane(Pnt(0, 0, Rminus), Vec(0, 0, -1)).bc('crack')
    pl2 = Plane(Pnt(0, 0, Rminus + delta), Vec(0, 0, 1))
    o_minus = (box - pl1).maxh(hsample)
    geometry.Add(o_minus.mat('ominus').maxh(hmax))
    geometry.Add((o_ext - pml).mat('pml'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - box).mat('oplus').maxh(hmax))
    geometry.Add((box * pl1 * pl2).mat('olayer').maxh(hmax))
    geometry.CloseSurfaces(pl1, pl2)
    return geometry


def sphere_geometry(Rminus, Rplus, Rext, Rpml, Rother, c, delta, hsample, hmax
    ):
    geometry = CSGeometry()
    o_ext = Sphere(Pnt(0, 0, 0), Rext).bc('outer')
    pml = Sphere(Pnt(0, 0, 0), Rpml)
    o_plus = Sphere(Pnt(0, 0, 0), Rplus).bc('interface')
    o_minus = Sphere(Pnt(0, 0, 0), Rminus).maxh(hmax).mat('ominus').maxh(
        hsample)
    other = Sphere(Pnt(0, c, 0), Rother)
    withCrack = o_minus * other
    withoutCrack = o_minus - other
    geometry.Add((o_ext - pml).mat('air'))
    geometry.Add((pml - o_plus).mat('air'))
    geometry.Add((o_plus - o_minus - other).mat('oplus'))
    geometry.Add((other - o_minus).mat('olayer'))
    geometry.Add(withCrack, bcmod=[(o_minus, 'crack')])
    geometry.Add(withoutCrack, bcmod=[(o_minus, 'nocrack')])
    return geometry

from ngsolve import *
from netgen.csg import *
from ngsolve.internal import *


def ATCerror_brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    box = OrthoBrick(Pnt(-Rminus,-Rminus,-Rminus),Pnt(Rminus,Rminus,Rminus+delta))
    pl1 = Plane(Pnt(0,0,Rminus),Vec(0,0,-1)).bc("crack")
    pl2 = Plane(Pnt(0,0,Rminus+delta),Vec(0,0,1))#.bc("top")
    o_minus = (box - pl1)
    geometry.Add ((box - pl1).mat("ominus"),bcmod=[(o_minus,"nocrack")])    
    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-box).mat("oplus").maxh(hmax))
    geometry.Add ((box * pl1 * pl2).mat("olayer").maxh(hmax),bcmod=[(pl1,"crack"),(box,"sides"),(pl2,"top")])
    #slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2)#,slices)
    
    return geometry

def ATCerror_halfsphere_geometry(Rminus, Rplus, Rext, Rpml, delta, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus+.2).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    circle = Sphere(Pnt(0,0,0),Rminus)
    small_cylinder = Cylinder(Pnt(0,0,1),Pnt(0,0,-1),Rminus*0.5)

    pl1 = Plane(Pnt(0,0,0),Vec(0,0,-1))
    pl2 = Plane(Pnt(0,0,delta),Vec(0,0,1))#.bc("top")

    o_minus_with_layer = (circle*pl2)
    o_minus = (circle - pl1).maxh(hmax)
    no_olayer = ((pl1 * pl2 * circle) - small_cylinder)
    yes_olayer =(pl1 * pl2 * small_cylinder)

    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-o_minus_with_layer).mat("oplus").maxh(hmax))
    geometry.Add ((o_minus).mat("ominus").maxh(hmax))    
    geometry.Add (yes_olayer.mat("olayer").maxh(hmax),bcmod=[(pl1,"crack")])
    geometry.Add (no_olayer.mat("oplus").maxh(hmax))
    slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2,slices)
    
    return geometry

def brick_geometry(Rminus, Rplus, Rext, Rpml, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")

    pml = Sphere(Pnt(0,0,0),Rpml)

    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    box = OrthoBrick(Pnt(-Rminus,-Rminus,-Rminus),Pnt(Rminus,Rminus,Rminus))
    pl1 = Plane(Pnt(0,0,Rminus),Vec(0,0,-1)).bc("crack")
    pl2 = Plane(Pnt(0,0,Rminus+delta),Vec(0,0,1))#.bc("top")
    o_minus = (box - pl1).maxh(hsample)

    geometry.Add (o_minus.mat("ominus").maxh(hmax))    
    geometry.Add ((o_ext - pml).mat("pml"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-box).mat("oplus").maxh(hmax))
    geometry.Add ((box * pl1 * pl2).mat("olayer").maxh(hmax))

    #slices = [2**(-i) for i in reversed(range(1,6))]
    geometry.CloseSurfaces(pl1,pl2)#,slices)
    
    return geometry

def sphere_geometry(Rminus, Rplus, Rext, Rpml, Rother, c, delta, hsample, hmax):
    geometry = CSGeometry()
    o_ext = (Sphere(Pnt(0,0,0), Rext)).bc("outer")
    pml = Sphere(Pnt(0,0,0),Rpml)
    o_plus = Sphere(Pnt(0,0,0), Rplus).bc("interface")

    #This is to define the two close surfaces for the thin layer:
    
    o_minus = (Sphere(Pnt(0,0,0), Rminus)).maxh(hmax).mat("ominus").maxh(hsample)
    other = Sphere(Pnt(0,c,0),Rother)
    withCrack = (o_minus * other)
    withoutCrack = (o_minus - other)

    geometry.Add ((o_ext - pml).mat("air"))
    geometry.Add ((pml-o_plus).mat("air"))
    geometry.Add ((o_plus-o_minus-other).mat("oplus"))
    geometry.Add ((other-o_minus).mat("olayer"))
    geometry.Add (withCrack,bcmod=[(o_minus,"crack")])
    geometry.Add (withoutCrack,bcmod=[(o_minus,"nocrack")])

    return geometry
