import os
import pip
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import sys
import math, random
import glob
from natsort import natsorted, ns
[docs]def force_install(package, versions=None):
"""install one package with versions """
if versions is not None:
pip.main(['install', package+'=='+versions])
else:
pip.main(['install', package])
[docs]def check_and_install():
try:
__import__('skgeom')
except ImportError:
print('...install... scikit-geometry')
cmd = 'conda install -c conda-forge scikit-geometry'
os.system(cmd)
check_and_install()
import skgeom as sg
from skgeom.draw import draw
"""
Generate random polygons with random BCs at random boundary locations for diffusion problem
"""
global TOTAL_COUNT
# https://stackoverflow.com/questions/8997099/algorithm-to-generate-random-2d-polygon
[docs]def generatePolygon( ctrX, ctrY, aveRadius, irregularity, spikeyness, numVerts ) :
'''Start with the centre of the polygon at ctrX, ctrY,
then creates the polygon by sampling points on a circle around the centre.
Randon noise is added by varying the angular spacing between sequential points,
and by varying the radial distance of each point from the centre.
Params:
ctrX, ctrY - coordinates of the "centre" of the polygon
aveRadius - in px, the average radius of this polygon, this roughly controls how large the polygon is, really only useful for order of magnitude.
irregularity - [0,1] indicating how much variance there is in the angular spacing of vertices. [0,1] will map to [0, 2pi/numberOfVerts]
spikeyness - [0,1] indicating how much variance there is in each vertex from the circle of radius aveRadius. [0,1] will map to [0, aveRadius]
numVerts - self-explanatory
Returns a list of vertices, in CCW order.
'''
irregularity = clip( irregularity, 0,1 ) * 2.0*math.pi / numVerts
spikeyness = clip( spikeyness, 0,1 ) * aveRadius
# generate n angle steps
angleSteps = []
lower = (2.*math.pi / numVerts) - irregularity
upper = (2.*math.pi / numVerts) + irregularity
sum = 0
for i in range(numVerts) :
tmp = random.uniform(lower, upper)
angleSteps.append( tmp )
sum = sum + tmp
# normalize the steps so that point 0 and point n+1 are the same
k = sum / (2.*math.pi)
for i in range(numVerts) :
angleSteps[i] = angleSteps[i] / k
# now generate the points
points = []
angle = random.uniform(0, 2*math.pi)
for i in range(numVerts) :
r_i = clip( random.gauss(aveRadius, spikeyness), 0, 2*aveRadius )
x = ctrX + r_i*math.cos(angle)
y = ctrY + r_i*math.sin(angle)
points.append( (int(x),int(y)) )
angle = angle + angleSteps[i]
return points
[docs]def clip(x, min, max) :
if( min > max ) : return x
elif( x < min ) : return min
elif( x > max ) : return max
else : return x
[docs]def create_polygon(points, pixel, shape_id=0, use_convex_hull=True):
sgPoints = []
for p0 in points:
sgPoints.append(sg.Point2(p0[0], p0[1]))
if use_convex_hull :
# the following make sure the inside has +1, and outside has -1
# but not work for L-shape, I guess.
chull_points = sg.convex_hull.graham_andrew(sgPoints)
poly = sg.Polygon(chull_points)
else:
# for L-shape, please define points counter-clock-wise
print("Please define the points counter-clock-wise to make sure inner has flag of +1 and outer has -1")
poly = sg.Polygon(sgPoints)
draw(sgPoints)
return poly
[docs]def generate_poly_points(num_points, pixel, use_new_random_polygon=False):
"""
generate points for polygons
"""
a = sg.random_polygon(num_points * 5, shape='circle', size=pixel/2)
a2 = sg.simplify(a, 0.2, "ratio", preserve_topology=False)
points = []
for v0 in list(a2.vertices):
x0 = int(float(v0.x())+pixel/2)
x0 = max(0,x0)
x0 = min(pixel-1, x0)
y0 = int(float(v0.y())+pixel/2)
y0 = max(0,y0)
y0 = min(pixel-1, y0)
p0 = [x0, y0]
points.append(p0)
if use_new_random_polygon:
while True:
verts = generatePolygon( ctrX=pixel/2, ctrY=pixel/2, aveRadius=pixel*0.45, irregularity=0.2, spikeyness=0.2, numVerts=5 )
verts = [[int(x[0]), int(x[1])] for x in verts]
list_of_verts = []
for x in verts:
list_of_verts.append(x[0])
list_of_verts.append(x[1])
# print(max(list_of_verts))
if max(list_of_verts) < pixel and min(list_of_verts) >= 0:
break
return verts
else:
return points
[docs]def distance(p1, p2, p0):
"""
line defined by p1 and p2
compute distance from p0 to the line
"""
d_list = []
_p1=np.array([float(p1.x()),float(p1.y())])
_p2=np.array([float(p2.x()),float(p2.y())])
_p0=np.array([float(p0.x()),float(p0.y())])
d_center = np.cross(_p2 - _p1, _p0 - _p1)/np.linalg.norm( _p2 - _p1 )
_d12 = np.linalg.norm( _p2 - _p1 )
_d10 = np.linalg.norm( _p1 - _p0 )
_d20 = np.linalg.norm( _p2 - _p0 )
d_list_max = -1.0
d_list_min = 1.0
if abs(d_center) < 0.5*np.sqrt(2):
_p0=np.array([float(p0.x())-0.5,float(p0.y())-0.5])
d_1 = np.cross(_p2 - _p1, _p0 - _p1)/np.linalg.norm( _p2 - _p1 )
_p0=np.array([float(p0.x())+0.5,float(p0.y())-0.5])
d_2 = np.cross(_p2 - _p1, _p0 - _p1)/np.linalg.norm( _p2 - _p1 )
_p0=np.array([float(p0.x())+0.5,float(p0.y())+0.5])
d_3 = np.cross(_p2 - _p1, _p0 - _p1)/np.linalg.norm( _p2 - _p1 )
_p0=np.array([float(p0.x())-0.5,float(p0.y())+0.5])
d_4 = np.cross(_p2 - _p1, _p0 - _p1)/np.linalg.norm( _p2 - _p1 )
d_list = [d_1, d_2, d_3, d_4]
d_list_max = max(d_list)
d_list_min = min(d_list)
# avoid potentially on the extension of the line
if (_d20 <= _d12 and _d10 <= _d12):
if d_list_max < 0: # four nodes are on one side of the line
return 999
elif d_list_min > 0: # four nodes are on one side of the line
return 999
else :
return d_center
else :
return 999
[docs]def plot_lines(pixel, color='b', linewidth=0.5):
for x in range(0, pixel):
plt.axhline(x+0.5, color=color, linestyle='-', linewidth=linewidth)
plt.axvline(x+0.5, color=color, linestyle='-', linewidth=linewidth)
[docs]def add_2edge_bc(bc_definition, count, total_edges):
for i in range(0, total_edges):
for j in range(i+2, i+2+total_edges-3):
# print(i, j%total_edges)
bc_definition[count] = {i:'c', j%total_edges:'h'}
count += 1
bc_definition[count] = {i:'c', j%total_edges:'c'}
count += 1
return count, bc_definition
# def add_3edge_bc(bc_definition, count, total_edges):
# for i in range(0, total_edges):
# for j in range(i+2, i+2+total_edges-5):
# for k in range(j+2, j+2+total_edges-5):
# # not loop back, for total_edges = 7+
# if (k+1)%total_edges != i:
# # not loop back, for total_edges = 8+
# if i != k%total_edges:
# # print(i, j%total_edges, k%total_edges)
# bc_definition[count] = {i:'c', j%total_edges:'c', k%total_edges:'c'}
# count += 1
# bc_definition[count] = {i:'c', j%total_edges:'c', k%total_edges:'h'}
# count += 1
# bc_definition[count] = {i:'c', j%total_edges:'h', k%total_edges:'h'}
# count += 1
# bc_definition[count] = {i:'c', j%total_edges:'h', k%total_edges:'c'}
# count += 1
# return count, bc_definition
[docs]def get_bc_definition(total_edges):
bc_definition = {}
if total_edges < 6:
bc_definition = {}
count = 0
count, bc_definition = add_2edge_bc(bc_definition, count, total_edges)
elif total_edges >= 6:
bc_definition = {}
count = 0
count, bc_definition = add_2edge_bc(bc_definition, count, total_edges)
# too much combinations
# count, bc_definition = add_3edge_bc(bc_definition, count, total_edges)
# print(bc_definition)
# exit(0)
else:
print("total_edges = ", total_edges, " is not implemented for get_bc_definition()")
exit(0)
print("total_bcs: ", count)
return bc_definition
[docs]def generate_bc_values(bc_size, bc_num, bc_range=[0.5, 1.0], bc_order_type="constant"):
order = 0
decimals = 2
if bc_order_type == "constant":
order = 0
elif bc_order_type == "linear":
order = 1
elif bc_order_type == "quadratic":
order = 2
elif bc_order_type == "sinusoidal":
order = 1
else :
print("bc_order_type = ", bc_order_type, " is not implemented in generate_bc_values()")
exit(0)
# np.random.seed(0) # cannot use fixed seed, as we loop over different edges, it could easily mess up things with the same value
values = np.random.uniform(bc_range[0], bc_range[1], bc_num * bc_size * (order+1) )
values = np.round(values, decimals)
values = values.reshape(bc_num, bc_size, (order+1))
list_of_bc_values = values
return list_of_bc_values
[docs]def assign_region_mask (img, poly, pixel):
for i0 in range(0, pixel):
for j0 in range(0, pixel):
d_list = []
for e0 in poly.edges:
d = distance(e0.point(0), e0.point(1), sg.Point2(i0, j0))
d_list.append(d)
# on one of the edges
if abs(min(d_list)) < 0.5*np.sqrt(2):
img[j0, i0, :] = 999
# not on the edges
else:
if poly.oriented_side(sg.Point2(i0, j0)) == sg.Sign.NEGATIVE:
img[j0, i0, :] = -1.0
else:
img[j0, i0, :] = 1.0
# assign value to the inner and outer region of the domain
img_c = img[:, :, 0:1]
img_h = img[:, :, 1:3]
img_c = np.where(img_c != 1, img_c, -2)
img_h = np.where(img_h != 1, img_h, 0)
img = np.concatenate((img_c, img_h), axis=2)
return img
[docs]def compute_bc_value(one_bc_value, bc_order_type, e0, i0, j0, bc_type):
normalization_factor = 2.0
if abs(float(e0.direction().dx())) >= abs(float(e0.direction().dy())):
is_x_dependent = True
else:
is_x_dependent = False
edge_vector = np.array([float(e0.direction().dx()), float(e0.direction().dy())])
edge_length = np.linalg.norm(edge_vector)
point_to_p0_vector = np.array([i0-float(e0.point(0).x()), j0-float(e0.point(0).y())])
projected_vector = np.dot(point_to_p0_vector, edge_vector)/np.dot(edge_vector, edge_vector) * edge_vector
point_to_p0_distance = np.linalg.norm(projected_vector)
norm = np.array([float(e0.direction().dy()), float(e0.direction().dx())])
norm = norm/np.linalg.norm(norm)
norm = abs(norm)
if bc_order_type == "constant":
new_bc_value = one_bc_value[0]
elif bc_order_type == "linear":
new_bc_value = one_bc_value[0] + (one_bc_value[1]-one_bc_value[0]) * point_to_p0_distance / edge_length
elif bc_order_type == "quadratic":
# min, next, max
one_bc_value.sort()
y = np.array(one_bc_value)
coor = np.array([0, edge_length, 0.5*edge_length])
new_bc_value = fit_polynomial_form(coor, y, point_to_p0_distance, order=2)
elif bc_order_type == "sinusoidal":
# https://en.wikipedia.org/wiki/Sine_wave
# y = A sin(w*t + phi), phi=0, A= max-min, w=2*pi, t = (coor_at - coor_1) / (coor_2 - coor_1)
# min, max
one_bc_value.sort()
t = point_to_p0_distance / edge_length
new_bc_value = one_bc_value[0] + (one_bc_value[1] - one_bc_value[0]) * np.sin(2*np.pi*t)
else :
# we can certainly add higher order terms
print("bc_order_type = ", bc_order_type, " is not implemented in compute_bc_value()")
exit(0)
if bc_type == 'c':
return new_bc_value
elif bc_type == 'h':
unscaled_bc_value = new_bc_value * normalization_factor - 0.5 * normalization_factor
unscaled_bc_value_split = unscaled_bc_value * norm
bc_value_split = unscaled_bc_value_split / normalization_factor + 0.5
bc_x = bc_value_split[0]
bc_y = bc_value_split[1]
return [bc_x, bc_y]
[docs]def apply_bcs_one_edge(img, list_of_edges, pixel, e0, _val, one_bc_value, bc_order_type):
"""
apply bcs on one edge
"""
if _val == 'c':
channel = 0
elif _val == 'h':
channel = 1
else:
print("bc val type = ", _val, ". Not sure what it is! In apply_bcs_one_edge()")
exit(0)
for i0 in range(0, pixel):
for j0 in range(0, pixel):
# only distance of pixels on the boundary will be calculated to save computational time
if img[j0,i0,channel] == 999:
d = distance(e0.point(0), e0.point(1), sg.Point2(i0, j0))
if abs(d) < 0.5 * np.sqrt(2):
# if d < 0.5 and d >= 0.0: # inside the domain
# if d > -0.5 and d <= 0.0: # outside the domain
# print(d)
if channel == 0:
img[j0, i0, channel] = compute_bc_value(one_bc_value, bc_order_type, e0, i0, j0, bc_type=_val)
if channel == 1:
# bc_x = compute_bc_value(one_bc_value, bc_order_type, e0, i0, j0, bc_type=_val)
[bc_x, bc_y] = compute_bc_value(one_bc_value, bc_order_type, e0, i0, j0, bc_type=_val)
img[j0, i0, channel] = bc_x
img[j0, i0, channel+1] = bc_y
[docs]def update_neumann_bc_value(list_of_bc_values, num_ind, neumann_ind, dirichlet_ind, bc_order_type):
"""
update Neumann BC values to make sure that the final results are not exceeding 1.0
"""
# print(list_of_bc_values[num_ind, :, :])
new_max_value = 1.0
for i0 in dirichlet_ind:
# print(list_of_bc_values[num_ind, i0, :])
new_max_value = min(new_max_value, 1.0-np.max(list_of_bc_values[num_ind, i0, :]) + 0.5)
# generate new bc values for Neumann BCs
new_range = [0.5, new_max_value]
new_bc_values = generate_bc_values(bc_size=len(dirichlet_ind), bc_num=1, bc_range=new_range, bc_order_type=bc_order_type)
_count = 0
for i0 in neumann_ind:
# print(i0, _count, new_bc_values)
list_of_bc_values[num_ind, i0, :] = new_bc_values[0, _count, :]
_count += 1
# print('new_bc_values:', new_bc_values, new_range, list_of_bc_values[num_ind, :, :])
[docs]def apply_bcs(poly, pixel, bc_num=5, bc_range=[0.5, 1.0], bc_order_type="constant", shape_id=0, selected_bcs=5):
"""
Apply constant, linear, quadratic, sinusoidal for the diffusion problem
"""
global TOTAL_COUNT
img_ref = np.zeros([pixel, pixel, 3])
img_ref = assign_region_mask(img_ref, poly, pixel)
#draw(poly, alpha=0.4)
total_edges = 0
list_of_edges = []
for e0 in poly.edges:
list_of_edges.append(e0)
# print(total_edges, e0)
total_edges += 1
bc_definition = get_bc_definition(total_edges)
selected_bc_num = min(selected_bcs, len(bc_definition))
print('total edges: ', total_edges, 'total bc: ', len(bc_definition), 'selected_bc_num: ', selected_bc_num)
# exit(0)
all_keys = np.arange(len(bc_definition))
np.random.shuffle(all_keys)
# np.random.shuffle(all_keys)
# print(all_keys)
data_count = 0
# for key, val in bc_definition.items():
for key in all_keys[0: selected_bc_num]:
val = bc_definition[key]
bc_size = len(val.items())
print ("working on BCs: ", key, val, bc_size)
list_of_bc_values = generate_bc_values(bc_size=bc_size, bc_num=bc_num, bc_range=bc_range, bc_order_type=bc_order_type)
# continue
# print(list_of_bc_values, np.shape(list_of_bc_values))
for num_ind in range(0, bc_num):
filename_prefix = "e"+str(total_edges) + "-s"+str(shape_id) + "-" + bc_order_type + "-" + str(data_count)
img = np.copy(img_ref)
print('****working on ****: total edge: ', total_edges, 'shape_id: ', shape_id, 'bc_order_type: ', bc_order_type, 'bc_value_num_id: ', num_ind, 'TOTAL remaining:', TOTAL_COUNT)
TOTAL_COUNT -= 1
with_neumann = False
# need to use list for displacement fields
neumann_ind = []
dirichlet_ind = []
it0 = 0
for _key, _val in val.items():
if _val == 'h':
with_neumann = True
neumann_ind.append(it0)
else:
dirichlet_ind.append(it0)
it0 += 1
if with_neumann: update_neumann_bc_value(list_of_bc_values, num_ind, neumann_ind, dirichlet_ind, bc_order_type)
it0 = 0
for _key, _val in val.items():
e0 = list_of_edges[_key]
one_bc_value = list_of_bc_values[num_ind, it0, :]
print('one_bc key, value: ', _key, _val, one_bc_value)
apply_bcs_one_edge(img, list_of_edges, pixel, e0, _val, one_bc_value, bc_order_type)
it0 += 1
# if with_neumann: print(list_of_bc_values[num_ind, :, :])
# create_cubit_mesh_dealii_bc(one_bc=val, one_bc_value=list_of_bc_values[num_ind, :, :], bc_order_type=bc_order_type, filename_prefix=filename_prefix)
# reset the free BC pixel to the correct mask value
img_c = img[:, :, 0:1]
img_h = img[:, :, 1:3]
img_c = np.where(img_c != 999, img_c, -2)
img_h = np.where(img_h != 999, img_h, 0)
img = np.concatenate((img_c, img_h), axis=2)
img = remove_extra_bc_pixel(img, pixel)
# not showing the pngs
if False :
plt.clf()
plt.subplot(131)
plt.imshow(img[:,:,0])
draw(poly, alpha=0.5)
plot_lines(pixel)
plt.subplot(132)
plt.imshow(img[:,:,1])
plot_lines(pixel)
draw(poly, alpha=0.5)
plt.subplot(133)
plt.imshow(img[:,:,2])
plot_lines(pixel)
draw(poly, alpha=0.5)
plt.savefig(filename_prefix+'.png')
plt.show()
# plt.close()
# exit(0)
feature = np.expand_dims(img, axis=0)
label = feature[:,:,:,0:1]
label = np.where(label != -2, label, 0.5)
feature_filename = "np-features" + "-" + filename_prefix + ".npy"
np.save(feature_filename, feature)
label_filename = feature_filename.replace("features", "labels")
np.save(label_filename, label)
data_count += 1
[docs]def vary_domain(all_num_points, pixel, bc_num=5, total_shapes=10, selected_bcs=5):
"""
randomly generate polygons and apply different BCs to them.
"""
# counter_clock_wise, inside +1, outside -1
for num_points in all_num_points:
print("num_points", num_points)
# how many shapes per group of shapes
for shape_id in range(0, total_shapes):
points = generate_poly_points(num_points, pixel)
poly = create_polygon(points, pixel, shape_id=shape_id)
# create_cubit_mesh(poly, pixel, shape_id=shape_id)
# draw(poly)
# plt.show()
apply_bcs(poly, pixel, bc_order_type="constant", bc_num=bc_num, shape_id=shape_id, selected_bcs=selected_bcs)
apply_bcs(poly, pixel, bc_order_type="linear", bc_num=bc_num, shape_id=shape_id, selected_bcs=selected_bcs)
apply_bcs(poly, pixel, bc_order_type="quadratic", bc_num=bc_num, shape_id=shape_id, selected_bcs=selected_bcs)
apply_bcs(poly, pixel, bc_order_type="sinusoidal", bc_num=bc_num, shape_id=shape_id, selected_bcs=selected_bcs)
[docs]def merge_one_folder(data_folder):
"""
Merge all the generated data into one file and delete individual files
"""
file_list = glob.glob(data_folder + '/np-features*.npy')
file_list = natsorted(file_list, alg=ns.IGNORECASE)
# print (file_list)
features = None
labels = None
count = 0
for f1 in file_list:
print('file: ', count, f1)
count += 1
one_feature = np.load(f1)
label_path = f1.replace('features', 'labels')
one_label = np.load(label_path)
print('file:', f1, 'label:', np.shape(one_label), 'feature:', np.shape(one_feature))
if (features is None):
features = np.copy(one_feature)
labels = np.copy(one_label)
else:
features = np.concatenate((features, one_feature), axis=0)
labels = np.concatenate((labels, one_label), axis=0)
feature_name = "np-features-all.npy"
label_name = "np-labels-all.npy"
np.save(feature_name, features)
np.save(label_name, labels)
print('feature shape: ', np.shape(features))
print('label shape: ', np.shape(labels))
cmd = 'rm np-*constant*.npy np-*linear*.npy np-*quadratic*.npy np-*sinusoidal*.npy -rf'
os.system(cmd)
if __name__ == "__main__":
# example run
Pixel = 32
BC_num = 5
Total_shape = 5
Selected_BCs = 4
Types_BCs = 4 # 4 types of BCs: constant, linear, quadratic, sinusoidal
TOTAL_COUNT = len(Edge_List) * Types_BCs * Total_shape * Selected_BCs * BC_num
print('total_count', TOTAL_COUNT)
vary_domain(all_num_points=Edge_List, pixel=Pixel, bc_num=BC_num, total_shapes=Total_shape, selected_bcs=Selected_BCs)
merge_one_folder('./')