kinnala / scikit-fem

Simple finite element assemblers
https://scikit-fem.readthedocs.io
BSD 3-Clause "New" or "Revised" License
508 stars 80 forks source link

enforcing Dirichlet conditions in ex21 raises IndexError or TypeError #609

Closed gdmcbain closed 3 years ago

gdmcbain commented 3 years ago

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

https://github.com/kinnala/scikit-fem/blob/56fe1564d39c6ae0ea94c98c4ace2f864e4db106/docs/examples/ex21.py#L88

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 commented 3 years ago

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 commented 3 years ago

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.

ex21-condense

gdmcbain commented 3 years ago

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

gdmcbain commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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

https://github.com/kinnala/scikit-fem/blob/56fe1564d39c6ae0ea94c98c4ace2f864e4db106/docs/examples/ex19.py#L81-L84

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