Error: Too many indicies for array - trying to output properties of surface facets

I/O
1

HI everyone,

I currently have a box mesh that I have created externally and then imported into Fenics, the files for which can be found below:
https://drive.google.com/drive/folders/1vVMnGd-HdyTaMUhDoLFt72c92kuewxFm?usp=sharing

From this box mesh, I have created a group called ‘surface’ in the .msh file which tags all facets that lie on the surface of the mesh. From this group, I want to output the the coordinates of the facets midpoint and the normals to that facet. In addition I would also like to output the value of my solution Temperature at this midpoint (perhaps via averaging the value of the solution at the dof points that border the facet?)

Currently I have this, which is designed to output the midpoints and normals of each facet on the surface

from dolfin import *
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
from collections import defaultdict


data = 'box'
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
hdf.read(mesh, '/mesh', False)
cells_f = MeshFunction('size_t', mesh, 3)
hdf.read(cells_f, '/cells')
facets_f = MeshFunction('size_t', mesh, 2)
hdf.read(facets_f, '/facets')

boundaries = MeshFunction('size_t', mesh, 'box_facet_region.xml') 


Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1] + 0.1*x[2]", degree=3), Space) # Temperature solution across mesh
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))
T_values = abs(Temperature.vector().get_local().reshape((-1,1))) # obtain the local temperature values at the correponding midpoint, unsure of how to complete???
V = FunctionSpace(mesh,'P', 1)

tag = 2  #Tag for 'surface' group in .msh file
surfacefacetindices = np.where(boundaries.array() == tag)[0]
print(surfacefacetindices)
print(np.array(facets(mesh)))
surfacefacets = np.array(facets(mesh))[surfacefacetindices]
for f in surfacefacets:
    print('midpoint: {}, normal: {}'.format(f.midpoint().array(), f.normal().array()))

However, I am currently getting an error with the definition of surfacefacets

File "outputtext_v2.py", line 60, in <module>
    surfacefacets = np.array(facets(mesh))[surfacefacetindices]
IndexError: too many indices for array

I am trying to pick out the corresponding surface facets from the list of surfacefacetindices but I am unsure how to resolve this problem. Does anyone know how I can do this and also output the corresponding values of Temperature at the midpoint of each facet?

Thanks in advance

2

If you inspect what you obtain by calling:

you will see that this is an array of an iterator, i.e.

In [1]: np.array(facets(mesh))                                                                                                                                                   
Out[1]: array(<dolfin.cpp.mesh.FacetIterator object at 0x7f98601d0d18>, dtype=object)

In [2]: np.array(facets(mesh)).shape                                                                                                                                             
Out[2]: ()

The fix is to wrap it as a list first:

surfacefacets = np.array(list(facets(mesh)))[surfacefacetindices]