Invalid value encountered in interpolate

dolfinx
1

Hi, I have a dirichlet boundary function in which a divide exists. The division is valid on the boundary but not in the whole domain. This is the warning:

RuntimeWarning: invalid value encountered in divide
  theta = np.arctan(x[1] / x[0])
2

Use the cells input to interpolate to limit the cells you are using interpolation on

3

Dear dokken,
My code is

x0 = 23.43211

y0 = 5.19

alpha = np.arctan(y0 / x0)

beta = np.arctan(x0 / y0)

def out_boundary(x):

    theta = np.arctan(x[1] / x[0])

    return 0.13395 - (theta - alpha) / (beta - alpha) * 0.2679

mesh_comm = MPI.COMM_WORLD

gdim = 3

gmsh_model_rank = 0

domain, cell_markers, facet_markers = gmshio.read_from_msh('pole-FEMMeshGmsh.msh', mesh_comm, rank=gmsh_model_rank, gdim=gdim)

#with XDMFFile(MPI.COMM_WORLD, "rfq_out.xdmf", "w") as xdmf:

# xdmf.write_mesh(domain)

# xdmf.write_meshtags(facet_markers)

V = FunctionSpace(domain, ("CG", 2))

u = TrialFunction(V)

v = TestFunction(V)

a = dot(grad(u), grad(v)) * dx

L = Constant(domain, ScalarType(0)) * v * dx

neg_pole_facets = facet_markers.indices[facet_markers.values == 27]

neg_pole_bdofs = locate_dofs_topological(V, 2, neg_pole_facets)

pos_pole_facets = facet_markers.indices[facet_markers.values == 26]

pos_pole_bdofs = locate_dofs_topological(V, 2, pos_pole_facets)

print(neg_pole_bdofs, pos_pole_bdofs)

out_facets = facet_markers.indices[facet_markers.values == 25]

out_bdofs = locate_dofs_topological(V, 2, out_facets)

u_O = Function(V)

u_O.interpolate(out_boundary)

u_O is the function I want to interpolate on the out_facets. Then how to modify the last line in my code? Thank you.

4

Please make a minimal reproducible example next time (using a built in mesh).
Consider the following:

import dolfinx

from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)

def f(x):
    return x[1]/x[0]


def left_facet_locator(x):
    return np.isclose(x[0], 1)


mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
left_facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, left_facet_locator)

cells = dolfinx.mesh.compute_incident_entities(mesh, left_facets, mesh.topology.dim-1, mesh.topology.dim)
u.interpolate(f, cells)

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u]) as vtx:
    vtx.write(0.0)

returning

image
image1448×1265 113 KB

1 Like
5

Thank you very much. :smiley: