MeshWorker in 1D

[Scott Miller]

Hi all,

I want to follow up on a question that was on the mailing list back in 2010 (I'm not exactly sure how to post a follow-up to the old list).  Has there been any work to make MeshWorker work for 1D problems?  

Thanks.

-Scott

[Scott Miller]

I want to follow up on a question that was on the mailing list back in 2010 (I'm not exactly sure how to post a follow-up to the old list).  Has there been any work to make MeshWorker work for 1D problems?  

Due to the lack of response, along with my own tests, I think the answers is that MeshWorker does not work for 1D problems.  

This begs the question, "What needs to be modified so that a (DG) code can be fully dimension independent in deal.ii?"  Is the real problem that FE values on 1D faces (points, dim=0) do not work?  What is the suggested approach for fixing this?  

Thanks,

-Scott

[Wolfgang Bangerth]

> This begs the question, "What needs to be modified so that a (DG) code can be
> fully dimension independent in deal.ii?" Is the real problem that FE values
> on 1D faces (points, dim=0) do not work? What is the suggested approach for
> fixing this?

FEFaceValues used to not run in 1d, but I think it does now. In any case, at
least it should.

Your question begs another question: what actually happens if you use
MeshWorker in 1d?

Best
W.


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

[Scott Miller]

FEFaceValues used to not run in 1d, but I think it does now. In any case, at
least it should.

Your question begs another question: what actually happens if you use
MeshWorker in 1d?

Well, the easy answer is the backtrace:

Program received signal EXC_BAD_ACCESS, Could not access memory.

Reason: KERN_INVALID_ADDRESS at address: 0x0000000000000000

boost::shared_ptr<dealii::FEValuesBase<1, 1> >::operator-> (this=0x0) at shared_ptr.hpp:424

424        BOOST_ASSERT(px != 0);

(gdb) bt

#0  boost::shared_ptr<dealii::FEValuesBase<1, 1> >::operator-> (this=0x0) at shared_ptr.hpp:424

#1  0x0000000100dbaaf1 in dealii::MeshWorker::IntegrationInfo<1, 1>::initialize_data (this=0x7fff5fbfe120, data=@0x7fff5fbfd6d0) at integration_info.templates.h:29

#2  0x0000000100014937 in dealii::MeshWorker::IntegrationInfoBox<1, 1>::initialize<dealii::Vector<double> > (this=0x7fff5fbfdc70, el=@0x7fff5fbfef30, mapping=@0x7fff5fbfea30, data=@0x7fff5fbfdbe0, block_info=0x0) at integration_info.h:1168

#3  0x000000010000fa61 in Parabolic::ParabolicProblem<1>::assemble_rhs (this=0x7fff5fbfea30, solution=@0x10c90aa00, residual=@0x10c90b9f0, assemble_jacobian=false) at /Users/smiller/Dropbox/ED_ComparisonStudy/ParabolicApprox/parabolic.cc:444

#4  0x000000010004a803 in boost::_mfi::mf3<void, Parabolic::ParabolicProblem<1>, dealii::Vector<double>&, dealii::Vector<double>&, bool>::operator() (this=0x7fff5fbff778, p=0x7fff5fbfea30, a1=@0x10c90aa00, a2=@0x10c90b9f0, a3=false) at mem_fn_template.hpp:393

#5  0x000000010004a748 in boost::_bi::list4<boost::_bi::value<Parabolic::ParabolicProblem<1>*>, boost::arg<1>, boost::arg<2>, boost::arg<3> >::operator()<boost::_mfi::mf3<void, Parabolic::ParabolicProblem<1>, dealii::Vector<double>&, dealii::Vector<double>&, bool>, boost::_bi::list3<dealii::Vector<double>&, dealii::Vector<double>&, bool&> > (this=0x7fff5fbff788, f=@0x7fff5fbff778, a=@0x7fff5fbfe418) at bind.hpp:457

#6  0x000000010004a652 in boost::_bi::bind_t<void, boost::_mfi::mf3<void, Parabolic::ParabolicProblem<1>, dealii::Vector<double>&, dealii::Vector<double>&, bool>, boost::_bi::list4<boost::_bi::value<Parabolic::ParabolicProblem<1>*>, boost::arg<1>, boost::arg<2>, boost::arg<3> > >::operator()<dealii::Vector<double>, dealii::Vector<double>, bool> (this=0x7fff5fbff778, a1=@0x10c90aa00, a2=@0x10c90b9f0, a3=@0x7fff5fbfe477) at bind_template.hpp:116

#7  0x000000010004a3b8 in boost::detail::function::void_function_obj_invoker3<boost::_bi::bind_t<void, boost::_mfi::mf3<void, Parabolic::ParabolicProblem<1>, dealii::Vector<double>&, dealii::Vector<double>&, bool>, boost::_bi::list4<boost::_bi::value<Parabolic::ParabolicProblem<1>*>, boost::arg<1>, boost::arg<2>, boost::arg<3> > >, void, dealii::Vector<double>&, dealii::Vector<double>&, bool>::invoke (function_obj_ptr=@0x7fff5fbff778, a0=@0x10c90aa00, a1=@0x10c90b9f0, a2=false) at function_template.hpp:153

#8  0x000000010000a68b in boost::function3<void, dealii::Vector<double>&, dealii::Vector<double>&, bool>::operator() (this=0x7fff5fbff770, a0=@0x10c90aa00, a1=@0x10c90b9f0, a2=false) at function_template.hpp:759

#9  0x00000001000093f6 in ExplicitTimeDiscretization<1>::compute_residual (this=0x7fff5fbff688, u_i_index=0, stage_index=0) at explicit_time_discretization.templates.h:208

#10 0x0000000100008382 in ExplicitTimeDiscretization<1>::advance_SSP_5_4 (this=0x7fff5fbff688) at explicit_time_discretization.templates.h:595

#11 0x0000000100006729 in ExplicitTimeDiscretization<1>::advance (this=0x7fff5fbff688) at explicit_time_discretization.templates.h:231

#12 0x000000010000594d in Parabolic::ParabolicProblem<1>::run (this=0x7fff5fbfea30) at /Users/smiller/Dropbox/ED_ComparisonStudy/ParabolicApprox/parabolic.cc:650

#13 0x00000001000041c8 in main () at /Users/smiller/Dropbox/ED_ComparisonStudy/ParabolicApprox/parabolic.cc:674

The interesting thing is the #2 line above:  The code around integration_info.h 1168 is 

     p = std_cxx1x::shared_ptr<VectorData<VECTOR, dim, sdim> >(new VectorData<VECTOR, dim, sdim> (boundary_selector));

     p->initialize(data);

     boundary_data = p;

     boundary.initialize_data(p);

-Scott 

[Scott Miller]

There are some instantiations in mesh worker/integration_info.h for dim=1 that I think can be removed.  They look like maybe they were implemented before some of the more recent updates for dim=1 support.

Once I comment them out, I get an error from "include/deal.II/grid/tria_accessor.templates.h" line 2054.  This is a constructor for TriaAccessor<0, 1, spacedim>, and the body of the function simply throws an "ExcInternalError()".    I do not understand exactly why this needs to occur.  I think getting past this issue will make meshworker function just fine with dim=1.

[Wolfgang Bangerth]

On 04/07/2013 02:45 PM, Scott Miller wrote:
> There are some instantiations in mesh worker/integration_info.h for dim=1 that
> I think can be removed. They look like maybe they were implemented before
> some of the more recent updates for dim=1 support.

Quite possible. Patch?

> Once I comment them out, I get an error from
> "include/deal.II/grid/tria_accessor.templates.h" line 2054. This is a
> constructor for TriaAccessor<0, 1, spacedim>, and the body of the function
> simply throws an "ExcInternalError()". I do not understand exactly why this
> needs to occur. I think getting past this issue will make meshworker function
> just fine with dim=1.

If you send a small testcase that shows this then I'd be happy to take a look.

Thanks

[Scott T. Miller]

There are some instantiations in mesh worker/integration_info.h for dim=1 that
I think can be removed.  They look like maybe they were implemented before
some of the more recent updates for dim=1 support.

Quite possible. Patch?

Sure; attached is a patch for meshworker/integration_info.h that simply removes

the dim=1 instantiations.  Warning:  I have thoroughly tested this, just with the example

below.

"include/deal.II/grid/tria_accessor.templates.h" line 2054.  This is a
constructor for TriaAccessor<0, 1, spacedim>, and the body of the function
simply throws an "ExcInternalError()".    I do not understand exactly why this
needs to occur.  I think getting past this issue will make meshworker function
just fine with dim=1.

If you send a small testcase that shows this then I'd be happy to take a look.

Attached is a simple DG method for a scalar advection equation.  It works fine

with dim=2, and seems to be ok with dim=1 if the TriaAccessor error is simply

commented out.  

With dim=1, you will get an error from line 251, constructing the object:

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

Thanks.

-Scott

[Wolfgang Bangerth]

Scott,

> Sure; attached is a patch for meshworker/integration_info.h that simply
> removes
> the dim=1 instantiations.

Thanks.

> Attached is a simple DG method for a scalar advection equation. It
> works fine
> with dim=2, and seems to be ok with dim=1 if the TriaAccessor error is
> simply
> commented out.

OK, I've fixed this error. Thanks for the testcase -- it's now also in
the testsuite!

Best
Wolfgang