Hello everyone,
I have a problem when defining the boundary conditions to solve a PDE in a 2D domain with dolfinx.
I have an array of data detect_value that contains the value of u, the solution of the PDE, in some points of the boundary, and an array detect_loc that contains the positions in which u takes those values.
For this reason I would like to impose the Dirichlet BCs as “If we are in detect_loc[i] take the value detect_value[i]”.
I tried a “smart” method and a “brute” method but both are not working.
u_bdry = lambda x: np.dot(detect_value, np.isclose(x, detect_loc).all(axis = 1))
u_bdry = lambda x: detect_value[0] * np.isclose(x, detect_loc[0,:])+ detect_value[1] * np.isclose(x, detect_loc[1,:]) + ...
detect_value[2] * np.isclose(x, detect_loc[2,:])+ detect_value[3] * np.isclose(x, detect_loc[3,:]) + ...
detect_value[4] * np.isclose(x, detect_loc[4,:])+ detect_value[5] * np.isclose(x, detect_loc[5,:]) + ...
detect_value[6] * np.isclose(x, detect_loc[6,:])+ detect_value[7] * np.isclose(x, detect_loc[7,:]) + ...
detect_value[8] * np.isclose(x, detect_loc[8,:])+ detect_value[9] * np.isclose(x, detect_loc[9,:]) + ...
detect_value[10] * np.isclose(x, detect_loc[10,:])+ detect_value[11] * np.isclose(x, detect_loc[11,:]) + ...
detect_value[12] * np.isclose(x, detect_loc[12,:])+ detect_value[13] * np.isclose(x, detect_loc[13,:]) + ...
detect_value[14] * np.isclose(x, detect_loc[14,:])+ detect_value[15] * np.isclose(x, detect_loc[15,:]) + ...
detect_value[16] * np.isclose(x, detect_loc[16,:])+ detect_value[17] * np.isclose(x, detect_loc[17,:]) + ...
detect_value[18] * np.isclose(x, detect_loc[18,:])+ detect_value[19] * np.isclose(x, detect_loc[19,:]) + ...
detect_value[20] * np.isclose(x, detect_loc[20,:])+ detect_value[21] * np.isclose(x, detect_loc[21,:]) + ...
detect_value[22] * np.isclose(x, detect_loc[22,:])+ detect_value[23] * np.isclose(x, detect_loc[23,:]) + ...
detect_value[24] * np.isclose(x, detect_loc[24,:])+ detect_value[25] * np.isclose(x, detect_loc[25,:]) + ...
detect_value[26] * np.isclose(x, detect_loc[26,:])+ detect_value[27] * np.isclose(x, detect_loc[27,:]) + ...
detect_value[28] * np.isclose(x, detect_loc[28,:])+ detect_value[29] * np.isclose(x, detect_loc[29,:]) + ...
detect_value[30] * np.isclose(x, detect_loc[30,:])+ detect_value[31] * np.isclose(x, detect_loc[31,:]) + ...
detect_value[32] * np.isclose(x, detect_loc[32,:])+ detect_value[33] * np.isclose(x, detect_loc[33,:]) + ...
detect_value[34] * np.isclose(x, detect_loc[34,:])+ detect_value[35] * np.isclose(x, detect_loc[35,:]) + ...
detect_value[36] * np.isclose(x, detect_loc[36,:])+ detect_value[37] * np.isclose(x, detect_loc[37,:]) + ...
detect_value[38] * np.isclose(x, detect_loc[38,:])+ detect_value[39] * np.isclose(x, detect_loc[39,:]) + ...
detect_value[40] * np.isclose(x, detect_loc[40,:])+ detect_value[41] * np.isclose(x, detect_loc[41,:]) + ...
detect_value[42] * np.isclose(x, detect_loc[42,:])+ detect_value[43] * np.isclose(x, detect_loc[43,:]) + ...
detect_value[44] * np.isclose(x, detect_loc[44,:])+ detect_value[45] * np.isclose(x, detect_loc[45,:]) + ...
detect_value[46] * np.isclose(x, detect_loc[46,:])+ detect_value[47] * np.isclose(x, detect_loc[47,:]) + ...
detect_value[48] * np.isclose(x, detect_loc[48,:])+ detect_value[49] * np.isclose(x, detect_loc[49,:]) + ...
detect_value[50] * np.isclose(x, detect_loc[50,:])+ detect_value[51] * np.isclose(x, detect_loc[51,:]) + ...
detect_value[52] * np.isclose(x, detect_loc[52,:])+ detect_value[53] * np.isclose(x, detect_loc[53,:]) + ...
detect_value[54] * np.isclose(x, detect_loc[54,:])+ detect_value[55] * np.isclose(x, detect_loc[55,:]) + ...
detect_value[56] * np.isclose(x, detect_loc[56,:])+ detect_value[57] * np.isclose(x, detect_loc[57,:]) + ...
detect_value[58] * np.isclose(x, detect_loc[58,:])+ detect_value[59] * np.isclose(x, detect_loc[59,:])
The first method gives me this error:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
File /usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py:318, in Function.interpolate.<locals>._interpolate(u, cells)
317 try:
--> 318 self._cpp_object.interpolate(u._cpp_object, cells)
319 except AttributeError:
AttributeError: 'function' object has no attribute '_cpp_object'
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
Input In [14], in <cell line: 36>()
35 # in our case rob is the non-homo dirichlet condition
36 if not flag_verification:
37 # zeros = lambda x: np.zeros(x[0].size)
---> 39 boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_bdry),
40 BoundaryCondition("Dirichlet", 2, u_bdry),
41 BoundaryCondition("Dirichlet", 3, u_bdry)]
Input In [14], in BoundaryCondition.__init__(self, type, marker, values)
5 if type == "Dirichlet":
6 u_D = fem.Function(V)
----> 7 u_D.interpolate(values)
8 facets = np.array(facet_tag.indices[facet_tag.values == marker])
9 #facets = ft.indices[ft.values == marker] doesn't work
File /usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py:336, in Function.interpolate(self, u, cells)
333 map = mesh.topology.index_map(mesh.topology.dim)
334 cells = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
--> 336 _interpolate(u, cells)
File /usr/lib/python3.9/functools.py:877, in singledispatch.<locals>.wrapper(*args, **kw)
873 if not args:
874 raise TypeError(f'{funcname} requires at least '
875 '1 positional argument')
--> 877 return dispatch(args[0].__class__)(*args, **kw)
File /usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py:320, in Function.interpolate.<locals>._interpolate(u, cells)
318 self._cpp_object.interpolate(u._cpp_object, cells)
319 except AttributeError:
--> 320 self._cpp_object.interpolate(u, cells)
Input In [9], in <lambda>(x)
---> 40 u_bdry = lambda x: np.dot(detect_value, np.isclose(x, detect_loc).all(axis = 1))
42 # u_bdry = lambda x: detect_value[0] * np.isclose(x, detect_loc[0,:])+ detect_value[1] * np.isclose(x, detect_loc[1,:]) + ...
43 # detect_value[2] * np.isclose(x, detect_loc[2,:])+ detect_value[3] * np.isclose(x, detect_loc[3,:]) + ...
44 # detect_value[4] * np.isclose(x, detect_loc[4,:])+ detect_value[5] * np.isclose(x, detect_loc[5,:]) + ...
(...)
72
73 ###
75 n = ufl.FacetNormal(domain)
File <__array_function__ internals>:5, in isclose(*args, **kwargs)
File /usr/local/lib/python3.9/dist-packages/numpy/core/numeric.py:2358, in isclose(a, b, rtol, atol, equal_nan)
2356 yfin = isfinite(y)
2357 if all(xfin) and all(yfin):
-> 2358 return within_tol(x, y, atol, rtol)
2359 else:
2360 finite = xfin & yfin
File /usr/local/lib/python3.9/dist-packages/numpy/core/numeric.py:2339, in isclose.<locals>.within_tol(x, y, atol, rtol)
2337 def within_tol(x, y, atol, rtol):
2338 with errstate(invalid='ignore'):
-> 2339 return less_equal(abs(x-y), atol + rtol * abs(y))
ValueError: operands could not be broadcast together with shapes (3,55020) (60,2)
While the second gives me this other error:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [9], in <cell line: 1>()
42 u_bdry = lambda x: detect_value[0] * np.isclose(x, detect_loc[0,:])+ detect_value[1] * np.isclose(x, detect_loc[1,:]) + ...
---> 43 detect_value[2] * np.isclose(x, detect_loc[2,:])+ detect_value[3] * np.isclose(x, detect_loc[3,:]) + ...
44 detect_value[4] * np.isclose(x, detect_loc[4,:])+ detect_value[5] * np.isclose(x, detect_loc[5,:]) + ...
45 detect_value[6] * np.isclose(x, detect_loc[6,:])+ detect_value[7] * np.isclose(x, detect_loc[7,:]) + ...
File <__array_function__ internals>:5, in isclose(*args, **kwargs)
File /usr/local/lib/python3.9/dist-packages/numpy/core/numeric.py:2355, in isclose(a, b, rtol, atol, equal_nan)
2352 dt = multiarray.result_type(y, 1.)
2353 y = asanyarray(y, dtype=dt)
-> 2355 xfin = isfinite(x)
2356 yfin = isfinite(y)
2357 if all(xfin) and all(yfin):
TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
I built my code based on this tutorial, if you have questions I can explain further.
Thanks in advance, I hope you have a great day!