Getting fancy with FEniCS 2019

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)

Is something like the following what you’re looking for?

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)

F = FunctionSpace(mesh,'DG',0)
f = Function(F)

# Indicator function for the circle.  Not quite the same as the MWE,
# but I didn't try to parse exactly what the logic was; it should be
# easy to adapt, though.
ell0 = 0.02
Px, Py = 0.5, 0.5
indicator = Expression("sqrt(pow(x[0]-Px,2)+pow(x[1]-Py,2))<dist? 1.0 : 0.0",
                       Px=Px,Py=Py,dist=4.0*ell0,degree=0)
f.interpolate(indicator)

# Note: Re-naming the function lets you re-use Paraview state files.
f.rename("f","f")
File_f = XDMFFile(outputFolder + "/f.xdmf")
File_f.write(f)
1 Like

Since i started recently to use fenics i am looking to alternative solutions and possibilities that it may offers (and ways to improve my coding skill as well), and that’s for sure an approach that i didn’t thought of. Thank you very much for the reply.