enforcing Dirichlet conditions in ex21 raises IndexError or TypeError

[gdmcbain]

While following up the discussion of enforce #591, in particular in relation to generalized eigenvalue problems, I noticed that it had been used in

scikit-fem/docs/examples/ex21.py

Line 88 in 56fe156

L, x = solve(*enforce(K, M, D=ib.find_dofs()['fixed']))

but then that running the example raises quite a lot of noise.

/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Exception in Tkinter callback
Traceback (most recent call last):
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/tkinter/__init__.py", line 1885, in __call__
    return self.func(*args)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/tkinter/__init__.py", line 806, in callit
    func(*args)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/_backend_tk.py", line 253, in idle_draw
    self.draw()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/backend_tkagg.py", line 9, in draw
    super(FigureCanvasTkAgg, self).draw()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py", line 407, in draw
    self.figure.draw(self.renderer)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/artist.py", line 41, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/figure.py", line 1863, in draw
    mimage._draw_list_compositing_images(
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/image.py", line 131, in _draw_list_compositing_images
    a.draw(renderer)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/artist.py", line 41, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 445, in draw
    sorted(self.collections,
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 446, in <lambda>
    key=lambda col: col.do_3d_projection(renderer),
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/art3d.py", line 669, in do_3d_projection
    self.update_scalarmappable()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/collections.py", line 855, in update_scalarmappable
    self._facecolors = self.to_rgba(self._A, self._alpha)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/cm.py", line 333, in to_rgba
    rgba = self.cmap(x, alpha=alpha, bytes=bytes)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/colors.py", line 610, in __call__
    rgba = lut[xa]
IndexError: arrays used as indices must be of integer (or boolean) type

Is the ComplexWarning due to the loss of symmetry mentioned in #591 on 2021-04-02? I note that in §8.4.3 ‘Keeping the essential nodes’ of Ern & Guermond (2004, cited in #591):

If the bilinear form a is symmetric, this technique has the apparent drawback of breaking the symmetry of the stiffness matrix. Actually, symmetry is not lost if an iterative method is used to solve the linear system.

[gdmcbain]

That was with matplotlib 3.3.2; on (pip install -U matplotlib) upgrading to 3.4.1, it's:

/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/skfem/visuals/matplotlib.py:39: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.
  ax = Axes3D(fig)
/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Exception in Tkinter callback
Traceback (most recent call last):
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/tkinter/__init__.py", line 1885, in __call__
    return self.func(*args)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/tkinter/__init__.py", line 806, in callit
    func(*args)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/_backend_tk.py", line 239, in idle_draw
    self.draw()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/backend_tkagg.py", line 9, in draw
    super().draw()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py", line 406, in draw
    self.figure.draw(self.renderer)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/artist.py", line 74, in draw_wrapper
    result = draw(artist, renderer, *args, **kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/artist.py", line 51, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/figure.py", line 2737, in draw
    mimage._draw_list_compositing_images(
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/image.py", line 132, in _draw_list_compositing_images
    a.draw(renderer)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/artist.py", line 51, in draw_wrapper
    return draw(artist, renderer, *args, **kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 485, in draw
    sorted(self.collections,
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 471, in do_3d_projection
    return artist.do_3d_projection()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/_api/deprecation.py", line 431, in wrapper
    return func(*inner_args, **inner_kwargs)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/mpl_toolkits/mplot3d/art3d.py", line 796, in do_3d_projection
    self.update_scalarmappable()
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/collections.py", line 926, in update_scalarmappable
    self._mapped_colors = self.to_rgba(self._A, self._alpha)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/cm.py", line 360, in to_rgba
    rgba = self.cmap(x, alpha=alpha, bytes=bytes)
  File "/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/matplotlib/colors.py", line 641, in __call__
    lut.take(xa, axis=0, mode='clip', out=rgba)
TypeError: Cannot cast array data from dtype('complex128') to dtype('int64') according to the rule 'safe'

and the shown figure is blank. (With 3.3.2 it was a uniform turquoise rather than the old graded plot that I remember.)

[gdmcbain]

Switching back to condense, I see that the ComplexWarning is still raised, along with some other matplotlib stuff

/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/skfem/utils.py:178: ComplexWarning: Casting complex values to real discards the imaginary part
  y[I] = X
/home/gmcbain/miniconda3/envs/skfem/lib/python3.9/site-packages/skfem/visuals/matplotlib.py:39: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6.  This is consistent with other Axes classes.
  ax = Axes3D(fig)

but the figure is as I remembered it.

[gdmcbain]

I think there are a few different things going on here.

[gdmcbain]

Reverting to the current enforce makes

pytest -k ex21 tests/test_examples.py 

fail (although not terribly significantly, though I'm wondering why it fails here and not in CI):

FAILED tests/test_examples.py::TestEx21::runTest - AssertionError: (50196.03548679974+0j) != 50194.51136114997 within 1 delta (1.5241256497683935 dif...

[gdmcbain]

The eigenvalues L look like enforce:

array([ 50196.0354868 +0.j, 251543.5420673 +0.j, 301209.62243364+0.j,
       598306.24657617+0.j, 729163.50693683+0.j])

and condense:

array([ 50194.94436948+0.j, 251542.76238838+0.j, 301210.3979456 +0.j,
       598305.23144065+0.j, 729165.63565666+0.j])

which I don't think is anything to worry about. I mean my previous fears about spurious eigenvalues don't seem to be relevant in the interesting part of the spectrum.

[gdmcbain]

I'm thinking that it should be possible to use enforce #591 #596 in ex19 to clean up a lot of stuff like

scikit-fem/docs/examples/ex19.py

Lines 81 to 84 in 56fe156

	u_new = np.zeros_like(u) # zero Dirichlet conditions
	_, b1 = condense(csr_matrix(A.shape), # ignore condensed matrix
	B @ u, u_new, D=boundary, expand=False)
	u_new[interior] = backsolve(b1)

in IVPs #531 but I haven't quite worked it out yet.