Skip to content

Commit

Permalink
merge upstream
Browse files Browse the repository at this point in the history
  • Loading branch information
@brighthe
brighthe committed Oct 17, 2024
2 parents b85af90 + b39eb42 commit 6369d06
Show file tree
Hide file tree
Showing 17 changed files with 102 additions and 1,839 deletions.
2 changes: 1 addition & 1 deletion app/FuelRodSim/HeatEquationData.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ def domain(self):
def duration(self):
return [0, 1]

def source(self,p,t):
def source(self, p, t):
return 0

def dirichlet(self, p, t):
Expand Down
2 changes: 1 addition & 1 deletion example/experiment/fem/maxwell_eq_n2fem_dirichlet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@
maxit = args.maxit
maxit = 3
dim = args.dim
dim = 3
#dim = 3
backend = args.backend
bm.set_backend(backend)
if dim == 2:
Expand Down
37 changes: 30 additions & 7 deletions example/opt/Levy_example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@
levy_steps_tensor = bm.array(levy_steps, device=device, dtype=bm.float32)
bround_steps_tensor = bm.array(bround_steps, device=device, dtype=bm.float32)

L[n * i : (i + 1) * n, :, :] = bm.cumsum(levy_steps_tensor, axis=0)
B[n * i : (i + 1) * n, :, :] = bm.cumsum(bround_steps_tensor, axis=0)
L[1 : 1 + n + 1, :, :] = bm.cumsum(levy_steps_tensor, axis=0)
B[1 : 1 + n + 1, :, :] = bm.cumsum(bround_steps_tensor, axis=0)

end_time = time.perf_counter()
running_time = end_time - start_time
Expand All @@ -38,15 +38,38 @@
L_cpu = L.cpu().numpy()
B_cpu = B.cpu().numpy()

for particle in range(0, 1):
L_x_vals = L_cpu[:, 0, particle]
L_y_vals = L_cpu[:, 1, particle]
B_x_vals = B_cpu[:, 0, particle]
B_y_vals = B_cpu[:, 1, particle]

plt.scatter(L_cpu[0, 0, 0], L_cpu[0, 0, 0], marker=MarkerStyle('*'), color='red', s=500)
plt.xticks([])
plt.yticks([])

for i in range(1, 16):
if i < 6:
step = 50
elif i < 11:
step = 150
else:
step = 300
L_x_vals = L_cpu[ : step * i, 0, 0]
L_y_vals = L_cpu[ : step * i, 1, 0]
B_x_vals = B_cpu[ : step * i, 0, 0]
B_y_vals = B_cpu[ : step * i, 1, 0]
plt.figure()
plt.plot(L_x_vals, L_y_vals, '-', color = 'blue')
plt.plot(B_x_vals, B_y_vals, '-', color = 'black')
plt.scatter(L_x_vals[0], L_y_vals[0], marker=MarkerStyle('*'), color='red', s=500)
plt.xticks([])
plt.yticks([])

particle = 0
L_x_vals = L_cpu[ : , 0, particle]
L_y_vals = L_cpu[ : , 1, particle]
B_x_vals = B_cpu[ : , 0, particle]
B_y_vals = B_cpu[ : , 1, particle]
plt.figure()
plt.plot(L_x_vals, L_y_vals, '-', color = 'blue')
plt.plot(B_x_vals, B_y_vals, '-', color = 'black')
plt.scatter(L_x_vals[0], L_y_vals[0], marker=MarkerStyle('*'), color='red', s=500)
plt.xticks([])
plt.yticks([])
plt.show()
2 changes: 1 addition & 1 deletion fealpy/fem/dirichlet_bc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,7 +230,7 @@ def apply_vector(self, vector: TensorLike, matrix: SparseTensor,
else:
if uh is None:
uh = bm.zeros_like(f)
uh = self.space.boundary_interpolate(gd, uh, self.threshold)
uh, _ = self.space.boundary_interpolate(gd, uh, self.threshold)

bd_idx = self.boundary_dof_index
f = f - A.matmul(uh[:])
Expand Down
57 changes: 29 additions & 28 deletions fealpy/functionspace/cm_conforming_fe_space.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,9 @@
_S = slice(None)

class CmConformingFESpace2d(FunctionSpace, Generic[_MT]):
def __init__(self, mesh: _MT, p: int, m: int, isCornerNode: bool):
def __init__(self, mesh: _MT, p: int, m: int, isCornerNode: bool, device=None):
assert(p>4*m)
self.device = device
self.mesh = mesh
self.p = p
self.m = m
Expand Down Expand Up @@ -72,7 +73,7 @@ def node_to_dof(self):
ndof = self.number_of_internal_dofs('node')

NN = mesh.number_of_nodes()
n2d = bm.arange(NN*ndof, dtype=self.itype).reshape(NN, ndof)
n2d = bm.arange(NN*ndof, dtype=self.itype, device=self.device).reshape(NN, ndof)
return n2d

def edge_to_internal_dof(self):
Expand All @@ -83,7 +84,7 @@ def edge_to_internal_dof(self):
mesh = self.mesh
NN = mesh.number_of_nodes()
NE = mesh.number_of_edges()
e2id = bm.arange(NN*ndof, NN*ndof+NE*eidof,dtype=self.itype).reshape(NE, eidof)
e2id = bm.arange(NN*ndof, NN*ndof+NE*eidof,dtype=self.itype, device=self.device).reshape(NE, eidof)
return e2id
def cell_to_internal_dof(self):
p = self.p
Expand All @@ -96,7 +97,7 @@ def cell_to_internal_dof(self):
NE = mesh.number_of_edges()
NC = mesh.number_of_cells()
c2id = bm.arange(NN*ndof + NE*eidof,
NN*ndof+NE*eidof+NC*cidof,dtype=self.itype).reshape(NC, cidof)
NN*ndof+NE*eidof+NC*cidof,dtype=self.itype, device=self.device).reshape(NC, cidof)
return c2id
def cell_to_dof(self):
p = self.p
Expand All @@ -114,11 +115,11 @@ def cell_to_dof(self):
cidof = self.number_of_internal_dofs('cell')
eidof = self.number_of_internal_dofs('edge')
ndof = self.number_of_internal_dofs('node')
c2d = bm.zeros((NC, ldof), dtype=self.itype)
c2d = bm.zeros((NC, ldof), dtype=self.itype, device=self.device)
e2id = self.edge_to_internal_dof()
n2d = self.node_to_dof()
c2d[:, ldof-cidof:] = bm.arange(NN*ndof + NE*eidof,
NN*ndof+NE*eidof+NC*cidof,dtype=self.itype).reshape(NC, cidof)
NN*ndof+NE*eidof+NC*cidof,dtype=self.itype, device=self.device).reshape(NC, cidof)
c2d[:, :ndof*3] = n2d[cell].reshape(NC,-1)
c2eSign = mesh.cell_to_face_sign()
for e in bm.arange(3):
Expand All @@ -136,7 +137,7 @@ def is_boundary_dof(self, threshold): #TODO:这个threshold 没有实现
p = self.p
m = self.m
gdof = self.number_of_global_dofs()
isBdDof = bm.zeros(gdof, dtype=bm.bool)
isBdDof = bm.zeros(gdof, dtype=bm.bool, device=self.device)
isBdEdge = self.mesh.boundary_face_flag() #TODO:ds 中没有edge
isBdNode = self.mesh.boundary_node_flag()
# 边界边
Expand All @@ -162,7 +163,7 @@ def coefficient_matrix(self):

NC = mesh.number_of_cells()
ldof = self.number_of_local_dofs('cell')
tem = bm.eye(ldof, dtype=self.ftype)
tem = bm.eye(ldof, dtype=self.ftype, device=self.device)
coeff = bm.tile(tem, (NC, 1, 1))

# 多重指标
Expand All @@ -184,7 +185,7 @@ def coefficient_matrix(self):
# 局部自由度 顶点
S02m = [S02m0, S02m1, S02m2]
for v in range(3):
flag = bm.ones(3, dtype=bm.bool)
flag = bm.ones(3, dtype=bm.bool, device=self.device)
flag[v] = False
for alpha in S02m[v]:
i = midx2num(alpha[None, :])
Expand Down Expand Up @@ -212,7 +213,7 @@ def coefficient_matrix(self):
N = bm.einsum('cfd, ced->cef', glambda, n)
S1m = [S1m0, S1m1, S1m2]
for de in range(3):
e = bm.ones(3, dtype=bm.bool)
e = bm.ones(3, dtype=bm.bool, device=self.device)
e[de] = False # e
for alpha in S1m[de]:
i = midx2num(alpha[None, :])
Expand All @@ -234,7 +235,7 @@ def coefficient_matrix(self):
# 全局自由度
node = mesh.entity('node')
cell = mesh.entity('cell')
Ndelat = bm.zeros((NC, 3, 2, 2), dtype=self.ftype)
Ndelat = bm.zeros((NC, 3, 2, 2), dtype=self.ftype, device=self.device)
t0 = node[cell[:, 2]] - node[cell[:, 1]]
t1 = node[cell[:, 0]] - node[cell[:, 2]]
t2 = node[cell[:, 1]] - node[cell[:, 0]]
Expand All @@ -245,25 +246,25 @@ def coefficient_matrix(self):
Ndelat[:, 2, 0] = t1
Ndelat[:, 2, 1] = -t0
ndof = self.number_of_internal_dofs('node')
coeff1 = bm.zeros((NC, 3*ndof, 3*ndof), dtype=self.ftype)
coeff1 = bm.zeros((NC, 3*ndof, 3*ndof), dtype=self.ftype, device=self.device)
symidx = [symmetry_index(2, r) for r in range(1, 2*m+1)]
# 边界自由度
isBdNode = mesh.boundary_node_flag()
isBdEdge = mesh.boundary_face_flag()
edge = mesh.entity('edge')[isBdEdge]
n = mesh.edge_unit_normal()[isBdEdge]
NN = isBdNode.sum()
coeff2 = bm.tile(bm.eye(3*ndof, dtype=self.ftype), (NC, 1, 1))
nodefram = bm.zeros((NN, 2, 2), dtype=self.ftype)
bdnidxmap = bm.zeros(len(isBdNode),dtype=self.itype)
bdnidxmap[isBdNode] = bm.arange(NN, dtype=self.itype)
coeff2 = bm.tile(bm.eye(3*ndof, dtype=self.ftype, device=self.device), (NC, 1, 1))
nodefram = bm.zeros((NN, 2, 2), dtype=self.ftype, device=self.device)
bdnidxmap = bm.zeros(len(isBdNode),dtype=self.itype, device=self.device)
bdnidxmap[isBdNode] = bm.arange(NN, dtype=self.itype, device=self.device)
nodefram[bdnidxmap[edge[:, 0]], 1] += 0.5*n
nodefram[bdnidxmap[edge[:, 1]], 1] += 0.5*n
nodefram[:, 0] = nodefram[:, 1]@bm.array([[0,1],[-1, 0]],
dtype=self.ftype)
dtype=self.ftype, device=self.device)
kk = 0
for v in range(3):
flag = bm.ones(3, dtype=bm.bool)
flag = bm.ones(3, dtype=bm.bool, device=self.device)
flag[v] = False #v^*
coeff1[:, ndof*v, ndof*v] = 1
cidx = isBdNode[cell[:, v]] & ~isCornerNode[cell[:, v]]
Expand Down Expand Up @@ -311,8 +312,8 @@ def boundary_interpolate(self, gD, uh, threshold=None):
coridx = bm.where(isCornerNode)[0]
isBdNode = mesh.boundary_node_flag()
isBdEdge = mesh.boundary_face_flag()
bdnidxmap = bm.zeros(len(isBdNode), dtype=bm.int32)
bdnidxmap[isBdNode] = bm.arange(isBdNode.sum(),dtype=bm.int32)
bdnidxmap = bm.zeros(len(isBdNode), dtype=bm.int32, device=self.device)
bdnidxmap[isBdNode] = bm.arange(isBdNode.sum(),dtype=bm.int32, device=self.device)
n2id = self.node_to_dof()[isBdNode]
e2id = self.edge_to_internal_dof()[isBdEdge]

Expand All @@ -322,12 +323,12 @@ def boundary_interpolate(self, gD, uh, threshold=None):
NE = len(edge)

n = mesh.edge_unit_normal()[isBdEdge]
nodefram = bm.zeros((NN, 2, 2), dtype=bm.float64) ## 注意！！！先 t 后 n
nodefram = bm.zeros((NN, 2, 2), dtype=bm.float64, device=self.device) ## 注意！！！先 t 后 n
nodefram[bdnidxmap[edge[:, 0]], 1] += 0.5*n
nodefram[bdnidxmap[edge[:, 1]], 1] += 0.5*n
nodefram[:, 0] = nodefram[:, 1]@bm.array([[0, 1], [-1,
0]],dtype=bm.float64)
nodefram[bdnidxmap[coridx]] = bm.tile(bm.eye(2,dtype=bm.float64), (len(coridx), 1, 1))
0]],dtype=bm.float64, device=self.device)
nodefram[bdnidxmap[coridx]] = bm.tile(bm.eye(2,dtype=bm.float64, device=self.device), (len(coridx), 1, 1))

# 顶点自由度
uh[n2id[:, 0]] = gD[0](node)
Expand Down Expand Up @@ -450,8 +451,8 @@ def interpolation(self, flist):
coridx = bm.where(isCornerNode)[0]
isBdNode = mesh.boundary_node_flag()
isBdEdge = mesh.boundary_face_flag()
bdnidxmap = bm.zeros(len(isBdNode), dtype=bm.int32)
bdnidxmap[isBdNode] = bm.arange(isBdNode.sum(),dtype=bm.int32)
bdnidxmap = bm.zeros(len(isBdNode), dtype=bm.int32, device=self.device)
bdnidxmap[isBdNode] = bm.arange(isBdNode.sum(),dtype=bm.int32, device=self.device)
n2id = self.node_to_dof()[isBdNode]
e2id = self.edge_to_internal_dof()[isBdEdge]

Expand All @@ -461,12 +462,12 @@ def interpolation(self, flist):
NE = len(edge)

n = mesh.edge_unit_normal()[isBdEdge]
nodefram = bm.zeros((NN, 2, 2), dtype=bm.float64) ## 注意！！！先 t 后 n
nodefram = bm.zeros((NN, 2, 2), dtype=bm.float64, device=self.device) ## 注意！！！先 t 后 n
nodefram[bdnidxmap[edge[:, 0]], 1] += 0.5*n
nodefram[bdnidxmap[edge[:, 1]], 1] += 0.5*n
nodefram[:, 0] = nodefram[:, 1]@bm.array([[0, 1], [-1,
0]],dtype=bm.float64)
nodefram[bdnidxmap[coridx]] = bm.tile(bm.eye(2,dtype=self.ftype), (len(coridx), 1, 1))
0]],dtype=bm.float64, device=self.device)
nodefram[bdnidxmap[coridx]] = bm.tile(bm.eye(2,dtype=self.ftype, device=self.device), (len(coridx), 1, 1))

k = 1;
for r in range(1, 2*m+1):
Expand Down
3 changes: 2 additions & 1 deletion fealpy/functionspace/functional.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,8 @@ def symmetry_index(d, r, dtype=None):
"""
@brief 将 d 维 r 阶张量拉长以后，其对称部分对应的索引和出现的次数
"""
symidx0 = bm.tensor(list(combinations_with_replacement(range(d), r)), dtype=dtype)
symidx0 = bm.tensor(list(combinations_with_replacement(range(d), r)),
dtype=dtype)
coe = bm.flip(d**bm.arange(r, dtype=dtype))
symidx = bm.einsum('ij,j->i', symidx0, coe)

Expand Down
78 changes: 0 additions & 78 deletions fealpy/iopt/ANT_TSP.py
View file

This file was deleted.

Loading

0 comments on commit 6369d06

Please sign in to comment.