hp and SolutionTransfer -- ExcInterpolationNotImplemented ?!

[Denis Davydov]

Hi all,

Apparently, a long discussion on this topic in 2015 here https://groups.google.com/forum/?fromgroups#!searchin/dealii/hp|sort:relevance/dealii/a7RipcMRI5s/r5pQOr9FSf0J is not enough for myself.

I came across an issue which I don't think I saw before.

--------------------------------------------------------
An error occurred in line <859> of file </Users/davydden/libs-sources/deal.ii/davydden/source/fe/fe.cc> in function
    virtual void dealii::FiniteElement<3, 3>::get_interpolation_matrix(const FiniteElement<dim, spacedim> &, FullMatrix<double> &) const [dim = 3, spacedim = 3]
The violated condition was:
    false
Additional information:
    (none)
--------------------------------------------------------


Aborting!
----------------------------------------------------
libc++abi.dylib: terminating with uncaught exception of type dealii::FiniteElement<3, 3>::ExcInterpolationNotImplemented:
--------------------------------------------------------
An error occurred in line <859> of file </Users/davydden/libs-sources/deal.ii/davydden/source/fe/fe.cc> in function
    virtual void dealii::FiniteElement<3, 3>::get_interpolation_matrix(const FiniteElement<dim, spacedim> &, FullMatrix<double> &) const [dim = 3, spacedim = 3]
The violated condition was:
    false
Additional information:
    (none)
--------------------------------------------------------

Backtraces is useless and I can't really check where it is called:

* thread #1: tid = 0x3ce73e, 0x00007fff98e0fdd6 libsystem_kernel.dylib`__pthread_kill + 10, queue = 'com.apple.main-thread', stop reason = signal SIGABRT
  * frame #0: 0x00007fff98e0fdd6 libsystem_kernel.dylib`__pthread_kill + 10
    frame #1: 0x00007fff98efb787 libsystem_pthread.dylib`pthread_kill + 90
    frame #2: 0x00007fff98d75420 libsystem_c.dylib`abort + 129
    frame #3: 0x00007fff978d285a libc++abi.dylib`abort_message + 266
    frame #4: 0x00007fff978f7c37 libc++abi.dylib`default_terminate_handler() + 243
    frame #5: 0x00007fff98401ba3 libobjc.A.dylib`_objc_terminate() + 124
    frame #6: 0x00007fff978f4d69 libc++abi.dylib`std::__terminate(void (*)()) + 8
    frame #7: 0x00007fff978f49f2 libc++abi.dylib`__cxa_rethrow + 99
    frame #8: 0x00000001000318a5 CO_arpack_04_hp.debug`main(argc=1, argv=0x00007fff5fbfee40) + 8261 at CO_arpack_04_hp.cc:91
    frame #9: 0x00007fff98ce1255 libdyld.dylib`start + 1

All that is happening for a simple p-only refinement of 16 elements in a mesh from FE_Q<3>(1) to FE_Q<3>(2).

I am quite certain that interpolation matrix should be there.

I will be trying to re-produce this in a small unit test, but meanwhile wanted to post here just in case i miss something obvious.

The calls are as follows:

/. std::vector<typename Triangulation<dim>::active_cell_iterator> p_cells = get_cells_to_be_p_refined(); // and clear h-refinement flag

/.SolutionTransfer<...> soltrans(...)

/.store a copy of solution vector

/.triangulation.prepare_coarsening_and_refinement ();

/.soltrans.prepare_for_coarsening_and_refinement (...);

/.triangulation.execute_coarsening_and_refinement (); 

and only now loop through stored p_cells and increment active_fe_index:

              typename hp::DoFHandler<dim>::active_cell_iterator

              dof_cell   (&triangulation,

                          p_cell[i]->level(),

                          p_cell[i]->index(),

                          &dof_handler);

              dof_cell->set_active_fe_index(dof_cell->active_fe_index()+1);

followed by setup_system() and soltrans.interpolate()

For debuging, i output cell->level(), cell->index() and cell->active_fe_index() and all looks reasonable.

number of elements, as expected, don't change and number of dofs increase.

The only think I can immediately think of is the fact that I store <active_cell_iterators> to Tria. 

I could have instead saved iterators of DoFHandler, but this should not make a difference.

Regards,

Denis.

[Timo Heister]

Can you manually call get_interpolation_matrix() for every FE that is
in your collection? Otherwise, a call stack would be really helpful...

On Wed, Feb 15, 2017 at 12:34 PM, Denis Davydov <davy...@gmail.com> wrote:
> Hi all,
>
> Apparently, a long discussion on this topic in 2015 here

> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_forum_-3Ffromgroups-23-21searchin_dealii_hp&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=hsfFxe8uG10gyO3GSkSZ_OizXppydztTY5cO7mVE5nI&s=XI9Ky1KSvCMoTrI2WWleOK6gb-f0zgWfRUAqaPIlqYA&e= |sort:relevance/dealii/a7RipcMRI5s/r5pQOr9FSf0J

> --
> The deal.II project is located at https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=hsfFxe8uG10gyO3GSkSZ_OizXppydztTY5cO7mVE5nI&s=SwW1cBzuQomk98aCyjDyn3E7ySe6-tOREnHLwczMBPU&e=
> For mailing list/forum options, see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_forum_dealii-3Fhl-3Den&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=hsfFxe8uG10gyO3GSkSZ_OizXppydztTY5cO7mVE5nI&s=mhD40ydQzsND5wbu7Qdeb0DPXF9UqiR-e2noic4_PkQ&e=
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_optout&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=hsfFxe8uG10gyO3GSkSZ_OizXppydztTY5cO7mVE5nI&s=ndPW6p7PbBQ2WNdRz6nRKv1vJpl_XrYZ9PQz_D4jlnE&e= .



--
Timo Heister
http://www.math.clemson.edu/~heister/

[Wolfgang Bangerth]

On 02/15/2017 10:34 AM, Denis Davydov wrote:
>
> Backtraces is useless and I can't really check where it is called:

That's because it's thrown by `AssertThrow`, which yields an exception
that propagates up the thread's stack until it kills the thread. That's
the place you see for the backtrace.

If you run it in gdb, say
catch throw
and it will set a breakpoint on all `throw` statements (including the
ones in `AssertThrow`). Then run your program and the debugger will stop
everytime something gets thrown. There may be places where something
gets thrown legitimately (and caught), and in those cases just
`continue` in gdb, until gdb stops in the place you care about in fe.cc,
line 859.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Denis Davydov]

Thanks Wolfgang and Timo for a swift reply.

By adding catches on assertion, i go into:

solution_transfer.cc:211:

fe[i].get_interpolation_matrix (fe[j], matrices(i,j));

with i=0 and j=1

Then I examined my FECollection closer I realized that I was using FE_Q<> wrapped into FE_Enriched<>

even when I did not intend to use PUM. That amounted 

to interpolation from  FESystem(FE_Q<dim>(1), 1, FE_Nothing<>, 1) to FESystem(FE_Q<dim>(2), 1, FE_Nothing<>, 1).

By using FE_Q<> directly, everything worked, as expected!

Rookie mistake, nothing is broken! ;-)

Regards,

Denis.

[Wolfgang Bangerth]

> By adding catches on assertion,

It took me a few years until I learned this trick (the `catch throw`
thing) but it turns out to be really useful.

Cheers