This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Evaluate outside the part of the mesh a given process has access to

0 votes

I wrote a function that calculates the "inverse Curl" of a function; that is, given the equation $\vec{B} = \nabla \times \vec{A}$ one would solve for $\vec{A}$ by calculating the "inverse Curl" of $\vec{B}$ which for the $x-$component of $\vec{A}$ gives
$$
A_x (x, y, z)= \int_0^1[tzB_y(tx,ty,tz)-tyB_z(tx,ty,tz)]\mathrm{d}t,
$$
and similarly for the $y-$ and $z-$components. The following code works for serial execution but fails in parallel because a given process will need to evaluate the function B at points that do not belong to the part of the mesh that it is "devoted to" (think how in order to calculate the value of the inverse Curl at any point, the value at the origin is needed to compute the integral above). I'll appreciate any ideas on how to get this working in parallel.
Best regards,
Miguel Valdez

import dolfin
import numpy

mesh = dolfin.UnitCubeMesh(10, 10, 10)
VV = dolfin.VectorFunctionSpace(mesh, 'CG', 1)

B = dolfin.interpolate(dolfin.Constant((0.0, 0.0, 1.0)), VV)
dolfin.plot(B, title="The field")

exact = dolfin.project(dolfin.Expression(("-0.5*x[1]", "0.5*x[0]", "0.0")), VV)
dolfin.plot(exact, title="The 'exact' potential")

def invCurl(f, x):
    dt = 1.0/1000.0
    grid = numpy.linspace(0.0, 1.0, 1001)
    u = numpy.zeros(3)
    result = numpy.zeros(3)
    for t in grid:
        f.eval(u, t*x)
        result[0] += (t*x[2]*u[1]-t*x[1]*u[2])*dt
        result[1] += (t*x[0]*u[2]-t*x[2]*u[0])*dt
        result[2] += (t*x[1]*u[0]-t*x[0]*u[1])*dt
    return result

coords = mesh.coordinates()

A = dolfin.Function(VV)
A_local = A.vector().get_local()
j = 0
for i in range(0, len(A_local)-2, 3):
    A_local[i], A_local[i+1], A_local[i+2] = invCurl(B, coords[j])
    j += 1

A.vector().set_local(A_local)
A.vector().apply("insert")

dolfin.plot(A, title="The calculated potential")
dolfin.interactive()
asked May 19, 2015 by miguelvaldez FEniCS Novice (140 points)

1 Answer

0 votes

This is an ongoing issue for dolfin - see the bitbucket.org pages... You can turn on allow_extrapolation but this will not give you meaningful results.

In parallel, you should be able to check whether a given Point is in the local part of the Mesh or not, using the BoundingBox collision detection algorithms, built into dolfin.
On that basis, you can decide whether or not to eval at that point, and possibly communicate the result to other processes (manually, using MPI).

answered May 20, 2015 by chris_richardson FEniCS Expert (31,740 points)
...