-
Notifications
You must be signed in to change notification settings - Fork 0
/
run-gen.py
100 lines (70 loc) · 1.9 KB
/
run-gen.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#!/usr/bin/env python
# coding: utf8
import numpy as np
from assemb import MultiTrace, checker
from time import time
import scipy.linalg as la
#meshname = "./geo/sphere-disjoint.msh"
#meshname = "./geo/ellipse-disjoint.msh"
meshname = "./geo/all.msh"
kRef = 0.5 * np.pi
# dd = [
# { 'name': '0',
# 'phys': 1,
# 'union': [-1, -2],
# },
# { 'name': 'A',
# 'phys': 4,
# 'union': 1,
# },
# { 'name': 'B',
# 'phys': 9,
# 'union': 2,
# },
# # { 'name': 'C',
# # 'phys': 3,
# # 'union': 3,
# # }
# ]
mtf = MultiTrace(kRef, meshname, dd)
At, X, J, iJ = mtf.tolinop()
shape = mtf.shape
A = 2.0 * At
A2 = A * iJ * A
Ce = 0.5 * J - At
Ci = 0.5 * J + At
Ce2 = Ce * iJ * Ce
Ci2 = Ci * iJ * Ci
x = np.random.rand(shape[0]) + 1j * np.random.rand(shape[0])
xr = np.random.rand(shape[0])
checker('A2 = J', A2, J, x)
checker('exterior Proj.', Ce2, Ce, x)
checker('interior Proj.', Ci2, Ci, x)
checker('error-Calderon with random [no-sense]', A, J, x)
#################################################
#################################################
def dir_data(x, normal, dom_ind, result):
result[0] = -np.exp( 1j * kRef * x[1])
def neu_data(x, normal, dom_ind, result):
result[0] = -1j * normal[1] * kRef * np.exp( 1j * kRef * x[1])
#################################################
#################################################
b = mtf.rhs(dir_data, neu_data)
M = A - X
print('')
print(mtf.shape, flush=True)
print('')
#################################################
#################################################
#################################################
iA = iJ * A * iJ
#################################################
Pjac = iA
E = mtf.upper()
Pgs = iA + iA * E * iA
CCi = iJ * (0.5 * J + At)
CCe = (0.5 * J - At)
B = J - X
BD = B * CCi + CCe
F = B * CCi
# Msigma = lambda sigma: (A - J) + sigma * (J - X)