Hei,
I would like to modify a vector field in numpy and i want to map it back so it can be used as a Function.
I figured out the way this should be done is by defining a UserExpression and projecting it on the mesh (see option 1 in the code). However, this can be quite slow, therefore it is recommended to do this in C++ (see option 2). This is kind of complicated and i have not figured out how to use other Functions in an expression like that.
Now I wondered if there is something like option3() or option4() (see code below)? These don’t work properly, it looks like the underlying petc array of a Function is stored in a different order than the mesh coordinates (see figure). I guess there is some kind of mapping that happens in compute_vertex_values().
Thank you
import math
import dolfin as fe
import pygmsh
import numpy as np
import matplotlib.pyplot as plt
def rectangle(lcar):
with pygmsh.geo.Geometry() as geom:
p = [geom.add_point([.20, .20], lcar),
geom.add_point([0.0, .0], lcar),
geom.add_point([2.2, .0], lcar),
geom.add_point([2.2, .41], lcar),
geom.add_point([0.0, .41], lcar)]
c = [geom.add_line(p[1], p[2]),
geom.add_line(p[2], p[3]),
geom.add_line(p[3], p[4]),
geom.add_line(p[4], p[1])]
ll1 = geom.add_curve_loop([c[0], c[1], c[2], c[3]])
s = [geom.add_plane_surface(ll1)]
geom.add_surface_loop(s)
msh = geom.generate_mesh()
mesh = create2Dmesh(msh)
return mesh
def create2Dmesh(msh):
"""
Helping function to create a 2D mesh for FEniCS from a gmsh.
important! Dont leave any unused points like the center of the circle in
the node list. FEniCS will crash!
"""
msh.prune_z_0()
nodes = msh.points[1:]
cells = msh.cells_dict["triangle"].astype(np.uintp)-1
mesh = fe.Mesh()
editor = fe.MeshEditor()
# point, interval, triangle, quadrilateral, hexahedron
editor.open(mesh, "triangle", 2, 2)
editor.init_vertices(len(nodes))
editor.init_cells(len(cells))
[editor.add_vertex(i, n) for i, n in enumerate(nodes)]
[editor.add_cell(i, n) for i, n in enumerate(cells)]
editor.close()
return mesh
def option1(p, mesh):
a = 4.
class P(fe.UserExpression):
def eval(self, values, x):
values[0] = math.sin(x[0]*a)
f0 = P(degree=1, element=V.ufl_element())
p.assign(fe.project(f0, mesh=mesh))
return
def option2(p, mesh):
ue_code = '''
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
class U : public dolfin::Expression
{
public:
double a;
void eval(Eigen::Ref<Eigen::VectorXd> values,
Eigen::Ref<const Eigen::VectorXd> x,
const ufc::cell& cell) const override
{
values[0] = sin(x[0]*a);
}
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<U, std::shared_ptr<U>, dolfin::Expression>
(m, "U")
.def(py::init<>())
.def_readwrite("a", &U::a);
}
'''
f0 = fe.CompiledExpression(fe.compile_cpp_code(ue_code).U(), a=4.0,
degree=1)
p.assign(fe.project(f0, mesh=mesh))
return
def option3(p, mesh):
a = 4.0
x, y = np.split(mesh.coordinates(), 2, 1)
p.vector().set_local(np.sin(x.ravel()*a))
return
def option4(p, mesh):
a = 4.0
x, y = np.split(mesh.coordinates(), 2, 1)
p.vector().vec().array = np.sin(x.ravel()*a)
return
if __name__ == "__main__":
mesh = rectangle(.1)
V = fe.FunctionSpace(mesh, 'P', 1)
for i, option in enumerate([option1, option2, option3, option4]):
p = fe.Function(V) # initialized with zeros
option(p, mesh)
ary1 = p.vector().vec().array
ary2 = p.compute_vertex_values(mesh)
print(i+1, np.allclose(ary1, ary2))
x, y = np.split(mesh.coordinates(), 2, 1)
x, y = x.ravel(), y.ravel()
fig, (ax1, ax2) = plt.subplots(2, sharex=True, sharey=True)
ax1.tricontourf(x.ravel(), y.ravel(), mesh.cells(),
ary1, levels=15)
ax2.tricontourf(x.ravel(), y.ravel(), mesh.cells(),
ary2, levels=15)
ax1.set_aspect("equal")
ax2.set_aspect("equal")
ax1.set_title("p.vector().vec().array")
ax2.set_title("p.compute_vertex_values(mesh)")
plt.suptitle("option {:.0f}".format(i+1))
plt.show()