Hi,
i wrote this code that assign the value to a costant scalar function based on the distance to a point. As you can see the code is by no means efficent nor elegant. How can i improve it?
Thanks
from dolfin import *
import numpy as np
import mshr as mr
outputFolder = "Test/"
#initialize mesh and function
L, H = 1., 2.
cell_size = 0.01
nel = int(L / cell_size)
geom = mr.Rectangle(Point(0, 0), Point(L, H))
mesh = mr.generate_mesh(geom, nel)
mesh_file = XDMFFile(outputFolder + "/mesh.xdmf")
mesh_file.write(mesh)
F = FunctionSpace(mesh,'DG',0)
f = Function(F)
f.vector()[:]=2
ell0 = 0.02
Px, Py = 0.5, 0.5
cells = mesh.cells()
list_vert = mesh.coordinates()
#element(s) containing the point
bbt = mesh.bounding_box_tree()
cells_coll = bbt.compute_collisions(Point(Px, Py))
cells_0 = []
mid_point_0 = []
for ii in cells_coll:
if Cell(mesh, ii).contains(Point(Px, Py)) is True:
mid_point_0.append(Cell(mesh, ii).midpoint())
cells_0.append(ii)
print(cells_0)
f.vector()[cells_0] = 0
#centroid of the element(s)
def Area(P):
return np.absolute((P[0]*(P[3]-P[5]) + P[2]*(P[5]-P[1]) + P[4]*(P[1]-P[3])) /2.)
if len(mid_point_0) == 1:
x0 = mid_point_0[0].x()
y0 = mid_point_0[0].y()
else:
x_i = []
y_i = []
A_i = []
for v in celle_0:
P = Cell(mesh, v).get_vertex_coordinates()
A_i.append(Area(P))
for jj, v in enumerate(mid_point_0):
x_i.append(mid_point_0[jj].x())
y_i.append(mid_point_0[jj].y())
x0 = np.inner(x_i,A_i)/sum(A_i)
y0 = np.inner(y_i,A_i)/sum(A_i)
G = [x0, y0]
#loop over the element and check critical distance
for ll,k in enumerate(cells):
P1 = list_vert[cells[ll,0]]
P2 = list_vert[cells[ll,1]]
P3 = list_vert[cells[ll,2]]
G_i = [(P1[0]+P2[0]+P3[0])/3 , (P1[1]+P2[1]+P3[1])/3]
d = np.sqrt((G_i[0] - G[0])**2 + (G_i[1]-G[1])**2)
if f.vector()[ll] != 0:
if d <= 4*ell0:
f.vector()[ll]=1
File_f = XDMFFile(outputFolder + "/f.xdmf")
File_f.write(f)