Hi, I think I didn't make myself clear in the last email. Apologize. I have a circle in my model ( center at (50,50), radius=25) as attached which represents the zero level set. Then I built a field with circular level sets which values vary smoothly from the center of the circle representing the zero level set, outwards along the radius ( the field has a zero value on the surface of circle). Then I calculated the curvature of the field everywhere in the mesh and output the values to a file called 'data', and also plot the contours of curvature. My expectation is the curvature should equal to (1/r) and I will have many concentric circles in my plot. Meanwhile, every point on the circle representing zero level set should have a curvature of around 0.04 ( since r = 25, 1/r=0.04) I have two problems now. 1. when I look at the 'data' file and checked some point which suppose to be at the surface e.g (50, 25), (62, 71), the results are like 0.0 or 0.7 (not 0.04 as expected). That's where I get confused. How can I reach my expectation? |
curvature1.pdf
Description: Adobe PDF document
2. The shape of my plot is quite wavy with many wiggles, not concentric rings as expected. Thanks very much. Regards, Yunbo |
from fipy import * import numpy as np import matplotlib.pyplot as plt
mesh = Grid2D(dx=1.0, dy=1.0, nx=100, ny=100)
var = DistanceVariable(name='Test',
mesh=mesh,
value=-1.0,
hasOld=1)
x, y = mesh.getCellCenters()
var.setValue(1.0, where = ((x-50.0)**2+(y-50.0)**2>625.0))
var.calcDistanceFunction()
curvature = (var.getFaceGrad()).getDivergence()
f = open('data','w')
count = 0
for item in curvature:
f.write(str(item) + '\t')
count += 1
if count%100 == 0:
f.write('\n')
curvature = curvature.reshape((100,100))
N = [i for i in range(100)]
fig = plt.figure(1, figsize=(8,8))
ax = fig.add_subplot(111)
ax.contour(N,N,curvature)
plt.savefig('curvature1.pdf')
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
