AndreaCodd
commented
2 years ago `import numpy as np from esys.escript import * from esys.oxley import Rectangle from esys.weipa import saveSilo
print() print(" test 1 Oxley MT") print()
core dimensions
cx = 40000. cy = 50000.
initial oxley square size
sqx = 40000. sqy = 40000.
buffer dimensions
bx = 100000.
by = 90000.
ba = 100000.
domain dimensions and ground level
Dx = cx+bx*2. Dy = cy+by+ba grLevel = 0.
dx = Dx/2. dy = -Dy+ba
core position
cL = -cx/2 #-cx/2. cR = cx/2. cT = grLevel cB = grLevel-cy
number of elements in each direction
nx = Dx/sqx ny = Dy/sqy
print("\n entire domain H = "+str(Dx)+" V = "+str(Dy)) print(" H in ",str(-dx),str(dx)) print(" V in ",str(dy),str(dy+Dy)) print('airground interface at z = ',grLevel)
domain = Rectangle(n0 = nx, n1 = ny, l0 = (-dx, dx), l1 = (dy,dy+Dy)) domain.dump('test1_1.silo') print("\ninitial squares hv = "+str(sqx)+" "+str(sqy)) print("made initial domain") print()
refine core
domain.setRefinementLevel(4) domain.refineRegion(x0=cL,x1=cR,y0=cB, y1=cT)
domain.dump('test1_2.silo')`
esys-escript
Domain is 240km240km and is initially setup so that the initial squares are 40km40km.
Horizontal coordinates are between -120km to 120km. Vertical coordinates are between -140km and 100km.
The core region is 40km by 50km (so x in [-20, 20], y in [-50, 0].
The initial mesh does not have a line at y = 0km or at -50km (the bottom of the core) nor vertical lines at -20km and 20km (the sides of the core.
If I refine the core region 4 levels the result is