Skip to content

Instantly share code, notes, and snippets.

@hemmer hemmer/mesh.py
Created Jan 18, 2016

Embed
What would you like to do?
generating a nonuniform mesh
#!/usr/bin/env python
from scipy import optimize
from scipy import integrate
import numpy as np
import sys
l, L = 0.02, 1
N = 128
mu, sigma = 0.5, 0.1
# density distrubtion function
def rho(x):
return 1.0 + np.tanh((x-(mu-sigma))/l) - np.tanh((x-(mu+sigma))/l)
# indefinite integral of the density distrubtion function
def rhoint(x):
return x + l*np.log(np.cosh((x - mu + sigma)/l)) - l*np.log(np.cosh((-x+mu+sigma)/l))
# first print density distribution function
for x in np.linspace(0, 1, N):
print x, rho(x)
print
# and find normalisation constant
theta = integrate.quad(rho, 0, L)[0]
print >> sys.stderr, "theta =", theta
def f(x, x0):
return rhoint(x[0]) - rhoint(x0) - theta/N
xs = np.zeros(N)
x0 = 0
for i in xrange(N):
x0 = optimize.fsolve(f, x0, args=(x0,))[0]
xs[i] = x0
print x0, np.sin(np.pi * x0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.