Skip to content

GitLab

Commit 60733ab7 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleanup of code, documentation and tests

parent 409f97e7
Hide whitespace changes
Inline Side-by-side
Showing with 367 additions and 947 deletions
CMakeLists.txt
View file @ 60733ab7
......@@ -26,7 +26,6 @@ add_subdirectory(src)
add_subdirectory(dune)
add_subdirectory(doc)
add_subdirectory(lib)
add_subdirectory(example)
add_subdirectory(cmake/modules)
# finalize the dune project, e.g. generating config.h etc.
......
doc/geometries/ellipsoid.py
View file @ 60733ab7
......@@ -2,72 +2,17 @@ from sympy import *
import numpy as np
import collections
x, y, z = symbols("x y z")
a, b, c = symbols("a b c")
# simplify entries of an array
def simplify_all(A):
if isinstance(A, collections.Iterable):
return list(map(lambda a: simplify_all(a), A))
else:
return simplify(A)
def add(A,B):
if isinstance(A, collections.Iterable):
return list(map(lambda a,b: add(a,b), A,B))
else:
return A + B
def negate(A):
if isinstance(A, collections.Iterable):
return list(map(lambda a: negate(a), A))
else:
return -A
x, y, z = symbols("x y z", real=True)
a, b, c = symbols("a b c", positive=True)
X = [x,y,z]
# normal vector
div = sqrt(b**4*c**4*x**2 + a**4*c**4*y**2 + a**4*b**4**z**2)
N = [b**2*c**2*x/div, a**2*c**2*y/div, a**2*b**2*z/div]
print("N = ")
print("return {")
for i in range(3):
print(" ",ccode(N[i]),",")
print("};")
print("")
# projection
P = [[ (1 if i==j else 0) - N[i]*N[j] for j in range(3)] for i in range(3)]
# parametrization
nrm_X = sqrt(x**2 + y**2 + z**2)
X0 = [x/nrm_X, y/nrm_X, z/nrm_X]
u = atan(X0[1]/X0[0])
v = acos(X0[2])
X1 = [a*cos(u)*sin(v), b*sin(u)*sin(v), c*cos(v)]
print("X1 = ")
print("return {")
for i in range(3):
print(" ",ccode(X1[i]),",")
print("};")
print("")
# jacobian of parametrization
J = [[diff(X1[i],X[j]) for i in range(3)] for j in range(3)]
print("J = ")
print("return {")
for i in range(3):
print(" {")
for j in range(3):
print(" ",ccode(simplify(J[i][j])),",")
print(" },")
print("};")
print("")
# implicit representation of the surface
phi = x**2/a**2 + y**2/b**2 + z**2/c**2 - 1
print("phi = ")
print("return ",ccode(simplify(phi)),";")
# gradient of phi
Jphi = [diff(phi,X[i]) for i in range(3)]
print("Jphi = ")
print("return {")
......@@ -75,143 +20,3 @@ for i in range(3):
print(" ",ccode(simplify(Jphi[i])),",")
print("};")
print("")
# P(F)_j
def project1(F):
return simplify_all([
np.sum([
P[i][j] * F[i]
for i in range(3)])
for j in range(3)])
# P(F)_kl
def project2(F):
return simplify_all([[
np.sum([np.sum([
P[i][k] * P[j][l] * F[i][j]
for j in range(3)]) for i in range(3)])
for l in range(3)] for k in range(3)])
# P(F)_lmn
def project3(F):
return simplify_all([[[
np.sum([np.sum([np.sum([
P[i][l] * P[j][m] * P[k][n] * F[i][j][k]
for k in range(3)]) for j in range(3)]) for i in range(3)])
for n in range(3)] for m in range(3)] for l in range(3)])
# P(F)_mnop
def project4(F):
return simplify_all([[[[
np.sum([np.sum([np.sum([np.sum([
P[i][m] * P[j][n] * P[k][o] * P[l][p] * F[i][j][k][l]
for l in range(3)]) for k in range(3)]) for j in range(3)]) for i in range(3)])
for p in range(3)] for o in range(3)] for n in range(3)] for m in range(3)])
# Euclidean gradient
def Grad0(f):
return simplify_all([
diff(f, X[i])
for i in range(3)])
# surface gradient (covariant derivative of scalars)
def grad0(f):
return project1(Grad0(f))
def Grad1(T):
return simplify_all([[
diff(T[i], X[j])
for j in range(3)] for i in range(3)])
# surface shape operator
#B = negate(project2(Grad1(N)))
# covariant derivative of vector field
def grad1(T):
return simplify_all(add(project2(Grad1(T)), [[
np.sum([
B[i][j] * T[k] * N[k]
for k in range(3)])
for j in range(3)] for i in range(3)]))
def Grad2(T):
return simplify_all([[[
diff(T[i][j], X[k])
for k in range(3)] for j in range(3)] for i in range(3)])
def grad2(T):
return simplify_all([[[
np.sum([np.sum([np.sum([
P[L1][I1] * P[L2][I2] * P[L3][K] * diff(T[L1][L2], X[L3])
for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
np.sum([np.sum([
B[K][I1] * P[J2][I2] * T[L][J2] * N[L]
for J2 in range(3)]) for L in range(3)]) +
np.sum([np.sum([
B[K][I2] * P[J1][I1] * T[J1][L] * N[L]
for J1 in range(3)]) for L in range(3)])
for K in range(3)] for I2 in range(3)] for I1 in range(3)])
def Grad3(T):
return simplify_all([[[[
diff(T[i][j][k], X[l])
for l in range(3)] for k in range(3)] for j in range(3)] for i in range(3)])
def grad3(T):
return simplify_all(add(project4(Grad3(T)), [[[[
np.sum([np.sum([np.sum([
B[K][I1] * P[J2][I2] * P[J3][I3] * T[L][J2][J3] * N[L]
for J2 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
np.sum([np.sum([np.sum([
B[K][I2] * P[J1][I1] * P[J3][I3] * T[J1][L][J3] * N[L]
for J1 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
np.sum([np.sum([np.sum([
B[K][I3] * P[J1][I1] * P[J2][I2] * T[J1][J2][L] * N[L]
for J1 in range(3)]) for J2 in range(3)]) for L in range(3)])
for K in range(3)] for I3 in range(3)] for I2 in range(3)] for I1 in range(3)]))
# normal-rotation of scalar field
def rot0(f):
return [diff(f*N[2],y) - diff(f*N[1],z),
diff(f*N[0],z) - diff(f*N[2],x),
diff(f*N[1],x) - diff(f*N[0],y)]
def Div1(F):
return diff(F[0],x) + diff(F[1],y) + diff(F[2],z)
def div1(F):
return Div1(project1(F))
# def div2(t):
# F = Matrix([div1(t.row(0).T), div1(t.row(1).T), div1(t.row(2).T)])
# return P*F
# div(T)_I1,I2
def div3(T):
return simplify_all([[
np.sum([np.sum([np.sum([np.sum([np.sum([
P[L1][I1] * P[L2][I2] * P[L3][K] * P[L4][K] * diff(T[L1][L2][L3],X[L4])
for K in range(3)]) for L4 in range(3)]) for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
np.sum([np.sum([
B[I3][I1] * T[L][I2][I3] * N[L] +
B[I3][I2] * T[I1][L][I3] * N[L] +
B[I3][I3] * T[I1][I2][L] * N[L]
for I3 in range(3)]) for L in range(3)])
for I2 in range(3)] for I1 in range(3)])
#p0 = simplify_all( rot0(x*y*z) ) # => vector
#print("p0 = ", p0)
# p1 = simplify_all( grad1(p0) ) # => 2-tensor
# print("p1 = ", p1)
# f = simplify_all( add(negate(div3(grad2(p1))), p1) )
# print("f = ", f)
#F = simplify( -div2(grad1(p0)) + p0 )
#print("F = ", F)
#print("F*N = ", simplify(F.dot(N)))
\ No newline at end of file
doc/geometries/torus.py
View file @ 60733ab7
......@@ -7,9 +7,13 @@ R, r = symbols("R r", positive=True)
X = [x,y,z]
# implicit representation of the surface
phi0 = sqrt(x**2 + y**2) - R
phi = phi0**2 + z**2 - r**2
print("phi = ")
print("return ",ccode(simplify(phi)),";")
# gradient of phi
Jphi = [diff(phi,X[i]) for i in range(3)]
print("Jphi = ")
print("return {")
......@@ -17,240 +21,3 @@ for i in range(3):
print(" ",ccode(simplify(Jphi[i])),",")
print("};")
print("")
Hphi = [[diff(Jphi[i],X[j]) for j in range(3)] for i in range(3)]
print("Hphi = ")
print("return {")
for i in range(3):
print(" {")
for j in range(3):
print(" ",ccode(simplify(Hphi[i][j])),",")
print(" },")
print("};")
print("")
JHJ = np.sum([np.sum([ Jphi[i] * Hphi[i][j] * Jphi[j] for j in range(3)]) for i in range(3)])
normJ2 = np.sum([Jphi[i]*Jphi[i] for i in range(3)])
trH = np.sum([Hphi[i][i] for i in range(3)])
mean_curvature = (JHJ - normJ2*trH) / (2*normJ2*sqrt(normJ2))
print("mean_curvature = ",ccode(simplify(mean_curvature)))
# simplify entries of an array
def simplify_all(A):
if isinstance(A, collections.Iterable):
return list(map(lambda a: simplify_all(a), A))
else:
return simplify(A) #.subs(phi0, r**2-z**2)
def add(A,B):
if isinstance(A, collections.Iterable):
return list(map(lambda a,b: add(a,b), A,B))
else:
return A + B
def negate(A):
if isinstance(A, collections.Iterable):
return list(map(lambda a: negate(a), A))
else:
return -A
# Euclidean gradient
def Grad0(f):
return simplify_all([
diff(f, X[i])
for i in range(3)])
# grad_phi = Grad0(phi)
# div = sqrt(grad_phi[0]**2 + grad_phi[1]**2 + grad_phi[2]**2)
# N = simplify_all([ grad_phi[i]/div for i in range(3) ])
sqrt_x2_y2 = sqrt(x**2 + y**2)
N = [ x - 2*x/sqrt_x2_y2, y - 2*y/sqrt_x2_y2, z ]
# print("N = ")
# print("return {")
# for i in range(3):
# print(" ",ccode(N[i]),",")
# print("};")
# print("")
#div = sqrt(b**4*c**4*x**2 + a**4*c**4*y**2 + a**4*b**4*z**2)
#N = [b**2*c**2*x/div, a**2*c**2*y/div, a**2*b**2*z/div]
# projection
P = [[ (1 if i==j else 0) - N[i]*N[j] for j in range(3)] for i in range(3)]
# print("P = ")
# print("return {")
# for i in range(3):
# print(" {")
# for j in range(3):
# print(" ",ccode(P[i][j]),",")
# print(" },")
# print("};")
# print("")
# P(F)_j
def project1(F):
return [
np.sum([
P[i][j] * F[i]
for i in range(3)])
for j in range(3)]
# P(F)_kl
def project2(F):
return [[
np.sum([np.sum([
P[i][k] * P[j][l] * F[i][j]
for j in range(3)]) for i in range(3)])
for l in range(3)] for k in range(3)]
# P(F)_lmn
def project3(F):
return [[[
np.sum([np.sum([np.sum([
P[i][l] * P[j][m] * P[k][n] * F[i][j][k]
for k in range(3)]) for j in range(3)]) for i in range(3)])
for n in range(3)] for m in range(3)] for l in range(3)]
# P(F)_mnop
def project4(F):
return [[[[
np.sum([np.sum([np.sum([np.sum([
P[i][m] * P[j][n] * P[k][o] * P[l][p] * F[i][j][k][l]
for l in range(3)]) for k in range(3)]) for j in range(3)]) for i in range(3)])
for p in range(3)] for o in range(3)] for n in range(3)] for m in range(3)]
# surface gradient (covariant derivative of scalars)
def grad0(f):
return project1(Grad0(f))
def Grad1(T):
return [[
diff(T[i], X[j])
for j in range(3)] for i in range(3)]
# surface shape operator
B = negate(project2(Grad1(N)))
# print("B = ")
# print("return {")
# for i in range(3):
# print(" {")
# for j in range(3):
# print(" ",ccode(B[i][j]),",")
# print(" },")
# print("};")
# print("")
# covariant derivative of vector field
def grad1(T):
return add(project2(Grad1(T)), [[
np.sum([
B[i][j] * T[k] * N[k]
for k in range(3)])
for j in range(3)] for i in range(3)])
def Grad2(T):
return [[[
diff(T[i][j], X[k])
for k in range(3)] for j in range(3)] for i in range(3)]
def grad2(T):
return [[[
np.sum([np.sum([np.sum([
P[L1][I1] * P[L2][I2] * P[L3][K] * diff(T[L1][L2], X[L3])
for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
np.sum([np.sum([
B[K][I1] * P[J2][I2] * T[L][J2] * N[L]
for J2 in range(3)]) for L in range(3)]) +
np.sum([np.sum([
B[K][I2] * P[J1][I1] * T[J1][L] * N[L]
for J1 in range(3)]) for L in range(3)])
for K in range(3)] for I2 in range(3)] for I1 in range(3)]
def Grad3(T):
return [[[[
diff(T[i][j][k], X[l])
for l in range(3)] for k in range(3)] for j in range(3)] for i in range(3)]
def grad3(T):
return add(project4(Grad3(T)), [[[[
np.sum([np.sum([np.sum([
B[K][I1] * P[J2][I2] * P[J3][I3] * T[L][J2][J3] * N[L]
for J2 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
np.sum([np.sum([np.sum([
B[K][I2] * P[J1][I1] * P[J3][I3] * T[J1][L][J3] * N[L]
for J1 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
np.sum([np.sum([np.sum([
B[K][I3] * P[J1][I1] * P[J2][I2] * T[J1][J2][L] * N[L]
for J1 in range(3)]) for J2 in range(3)]) for L in range(3)])
for K in range(3)] for I3 in range(3)] for I2 in range(3)] for I1 in range(3)])
# normal-rotation of scalar field
def rot0(f):
return [diff(f*N[2],y) - diff(f*N[1],z),
diff(f*N[0],z) - diff(f*N[2],x),
diff(f*N[1],x) - diff(f*N[0],y)]
def Div1(F):
return diff(F[0],x) + diff(F[1],y) + diff(F[2],z)
def div1(F):
return Div1(project1(F))
# def div2(t):
# F = Matrix([div1(t.row(0).T), div1(t.row(1).T), div1(t.row(2).T)])
# return P*F
# div(T)_I1,I2
def div3(T):
return [[
np.sum([np.sum([np.sum([np.sum([np.sum([
P[L1][I1] * P[L2][I2] * P[L3][K] * P[L4][K] * diff(T[L1][L2][L3],X[L4])
for K in range(3)]) for L4 in range(3)]) for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
np.sum([np.sum([
B[I3][I1] * T[L][I2][I3] * N[L] +
B[I3][I2] * T[I1][L][I3] * N[L] +
B[I3][I3] * T[I1][I2][L] * N[L]
for I3 in range(3)]) for L in range(3)])
for I2 in range(3)] for I1 in range(3)]
#p0 = simplify_all( rot0(z) ) # => vector # rot0(x*y*z)
#print("p0 = ", p0)
#p1 = simplify_all( grad1(p0) ) # => 2-tensor
# p1 = simplify_all( project2([[x,0,0],[0,y,0],[0,0,z]]) )
# print("p1 = ")
# print("return {")
# for i in range(3):
# print(" {")
# for j in range(3):
# print(" ",ccode(p1[i][j]),",")
# print(" },")
# print("};")
# print("")
# f = simplify_all( add(negate(div3(grad2(p1))), p1) )
# print("f = ")
# print("return {")
# for i in range(3):
# print(" {")
# for j in range(3):
# print(" ",ccode(f[i][j]),",")
# print(" },")
# print("};")
#F = simplify( -div2(grad1(p0)) + p0 )
#print("F = ", F)
#print("F*N = ", simplify(F.dot(N)))
\ No newline at end of file
doc/grids/ellipsoid3.geo 0 → 100644
View file @ 60733ab7
Point(1) = {0, 0, 0, 0.75};
Point(2) = {1, 0, 0, 0.75};
Point(3) = {-1, 0, 0, 0.75};
Point(4) = {0, 0.75, 0, 0.75};
Point(5) = {0, -0.75, 0, 0.75};
Point(6) = {0, 0, 1.25, 0.75};
Point(7) = {0, 0, -1.25, 0.75};
Ellipse(1) = {2, 1, 2, 4};
Ellipse(2) = {4, 1, 3, 3};
Ellipse(3) = {3, 1, 3, 5};
Ellipse(4) = {5, 1, 2, 2};
Ellipse(5) = {2, 1, 6, 6};
Ellipse(6) = {6, 1, 6, 3};
Ellipse(7) = {3, 1, 7, 7};
Ellipse(8) = {7, 1, 7, 2};
Ellipse(9) = {4, 1, 6, 6};
Ellipse(10) = {6, 1, 6, 5};
Ellipse(11) = {5, 1, 7, 7};
Ellipse(12) = {7, 1, 7, 4};
Line Loop(1) = {2, 7, 12};
Surface(1) = {-1};
Line Loop(2) = {2, -6, -9};
Surface(2) = {2};
Line Loop(3) = {3, -10, 6};
Surface(3) = {3};
Line Loop(4) = {3, 11, -7};
Surface(4) = {-4};
Line Loop(5) = {4, -8, -11};
Surface(5) = {-5};
Line Loop(6) = {4, 5, 10};
Surface(6) = {6};
Line Loop(7) = {1, 9, -5};
Surface(7) = {7};
Line Loop(8) = {1, -12, 8};
Surface(8) = {-8};
doc/grids/sphere.geo 0 → 100644
View file @ 60733ab7
Point(1) = {0, 0, 0, 0.75};
Point(2) = {1, 0, 0, 0.75};
Point(3) = {-1, 0, 0, 0.75};
Point(4) = {0, 1, 0, 0.75};
Point(5) = {0, -1, 0, 0.75};
Point(6) = {0, 0, 1, 0.75};
Point(7) = {0, 0, -1, 0.75};
Ellipse(1) = {2, 1, 2, 4};
Ellipse(2) = {4, 1, 3, 3};
Ellipse(3) = {3, 1, 3, 5};
Ellipse(4) = {5, 1, 2, 2};
Ellipse(5) = {2, 1, 6, 6};
Ellipse(6) = {6, 1, 6, 3};
Ellipse(7) = {3, 1, 7, 7};
Ellipse(8) = {7, 1, 7, 2};
Ellipse(9) = {4, 1, 6, 6};
Ellipse(10) = {6, 1, 6, 5};
Ellipse(11) = {5, 1, 7, 7};
Ellipse(12) = {7, 1, 7, 4};
Line Loop(1) = {2, 7, 12};
Surface(1) = {-1};
Line Loop(2) = {2, -6, -9};
Surface(2) = {2};
Line Loop(3) = {3, -10, 6};
Surface(3) = {3};
Line Loop(4) = {3, 11, -7};
Surface(4) = {-4};
Line Loop(5) = {4, -8, -11};
Surface(5) = {-5};
Line Loop(6) = {4, 5, 10};
Surface(6) = {6};
Line Loop(7) = {1, 9, -5};
Surface(7) = {7};
Line Loop(8) = {1, -12, 8};
Surface(8) = {-8};
doc/grids/torus2.bt3 0 → 100644
View file @ 60733ab7
dim: 3
boxes: 3 3 1
box-size: 1.8 1.8 1.8
offset: -2.7 -2.7 -0.9
matrices:
z = 0:
1 1 1
1 0 1
1 1 1