Fix bug in MappingFEField #6296

luca-heltai · 2018-04-23T15:41:13Z

~~I'll commit here a bug fix that closes #6294.~~

~~For the moment it's only a failing test.~~

bangerth · 2018-04-23T16:28:54Z

include/deal.II/grid/grid_tools.h


  /**
   * Return the diameter of the largest active cell of a triangulation.
   */
  template <int dim, int spacedim>
  double
-  maximal_cell_diameter (const Triangulation<dim, spacedim> &triangulation);
+  maximal_cell_diameter (const Triangulation<dim, spacedim> &triangulation,
+                         const Mapping<dim,spacedim> &mapping = (StaticMappingQ1<dim,spacedim>::mapping));


I think you should document for these functions that the diameter is computed only by looking at the distances between vertices. For very deformed cells, this is not the diameter of the area/volume occupied by this cell.

I have not changed the original implementation. That's exactly what the function was doing before (i.e., that's what cell->diameter does). I'll document this as well.

BTW, the comments you made should really be on #6295. This branch is based on that one. I'll rebase afterwards.

davydden · 2018-04-23T18:30:54Z

source/grid/grid_tools.cc

+  namespace
+  {
+    template <int dim, int spacedim>
+    double diameter(const typename Triangulation<dim, spacedim>::cell_iterator &cell,


is this exactly what cell->diameter() does, ie. look at vertices? If so, can we avoid somehow code duplication and join the two?

davydden · 2018-04-23T18:31:40Z

source/grid/grid_tools.cc

+                                    (vertices[3]-vertices[4]).norm()) );
+        default:
+          Assert (false, ExcNotImplemented());
+          return -1e10;


std::numeric_limits<double>::signaling_NaN()

davydden · 2018-04-23T18:32:40Z

tests/grid/grid_tools_min_max_diameter_01.cc

+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+


add a quick note what this test does. Everyone will forget in half a year 😄

same for the other test.

davydden · 2018-04-24T12:31:50Z

include/deal.II/fe/mapping_fe_field.h

+  virtual
+  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const override;
+


we usually keep just one space for headers

I kept was was already there. Do you want me to fix all the rest of the header as well?

yeah, i saw that the header is different to whitespaces elsewhere. Up to you, i would probably just fix it to be one space everywhere.

davydden · 2018-04-24T12:33:05Z

source/fe/mapping_fe_field.cc

+namespace
+{
+  template<int dim>
+  Quadrature<dim> get_vertex_quadrature()


add a quick doxygen documentation

This is unnecessary. Its in the source and an anonymous namespace - its never exposed to the user.

I never said that it's for users 😃

Documentation of internal functions are also good to have so that we don't need to guess what's the purpose. Otherwise we might as well get rid of all in-code comments, no users will ever read them.

Otherwise we might as well get rid of all in-code comments, no users will ever read them.

Sigh. This was not particularly constructive.

Documentation of internal functions are also good to have so that we don't need to guess what's the purpose.

So, you recommend an inline comment to explain the purpose of this function, rather than making generated documentation?

This was not particularly constructive.

the point was that we do document internal functions/algorithms.

you recommend an inline comment to explain the purpose of this function

Does not matter which format it is (as it's never rendered) but a short one sentence explaining the function can never harm.

davydden · 2018-04-24T12:34:40Z

source/fe/mapping_fe_field.cc

+  typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell,
+                                                            euler_dof_handler);
+
+  Assert (dof_cell->active() == true, ExcInactiveCell());


you don't need == true

I prefer to be explicit.

Reading (dof_cell->active()) is worse than (dof_cell->active == true) which is what we do everywhere else...

Also because active is not, unfortunately, a good member name. is_active would have been more informative, and I'd have agreed with you then...

(in other words, active does not suggest what the return type may be, so == true makes this more informative)

fair enough.

p.s. is_active() would indeed be better...

davydden · 2018-04-24T13:35:16Z

/run-tests

davydden · 2018-04-24T13:35:40Z

Let's see if @bangerth has further comments.

bangerth · 2018-04-25T03:44:14Z

There is a problem:

161/3560 Test #163: mappings/mapping_fe_field_02.debug ...................................***Failed   17.64 sec
Test mappings/mapping_fe_field_02.debug: RUN
===============================   OUTPUT BEGIN  ===============================
[ 50%] Generating mapping_fe_field_02.debug/interrupt_guard.cc
Scanning dependencies of target mapping_fe_field_02.debug
[ 50%] Building CXX object CMakeFiles/mapping_fe_field_02.debug.dir/mapping_fe_field_02.cc.o
Building CXX object CMakeFiles/mapping_fe_field_02.debug.dir/mapping_fe_field_02.debug/interrupt_guard.cc.o
Linking CXX executable mapping_fe_field_02.debug/mapping_fe_field_02.debug
[100%] Built target mapping_fe_field_02.debug
Scanning dependencies of target mapping_fe_field_02.debug.diff
[100%] Generating mapping_fe_field_02.debug/output
mappings/mapping_fe_field_02.debug: BUILD successful.
mappings/mapping_fe_field_02.debug: RUN failed. ------ Return code 134
mappings/mapping_fe_field_02.debug: RUN failed. ------ Result: /home/bob/build-gcc/tests/mappings/mapping_fe_field_02.debug/failing_output
mappings/mapping_fe_field_02.debug: RUN failed. ------ Partial output:
JobId e68cde002913 Tue Apr 24 20:26:11 2018
DEAL::dim, spacedim: 1, 1
DEAL::cells: 1, dofs: 2

mappings/mapping_fe_field_02.debug: RUN failed. ------ Additional output on stdout/stderr:


--------------------------------------------------------
An error occurred in line <104> of file </home/bob/source/source/base/subscriptor.cc> in function
    void dealii::Subscriptor::check_no_subscribers() const
The violated condition was: 
    counter == 0
Additional information: 
    (none)

bangerth · 2018-04-25T03:48:43Z

source/fe/mapping_fe_field.cc

+    std::vector<Point<dim> > points(GeometryInfo<dim>::vertices_per_cell);
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      points[i] = GeometryInfo<dim>::unit_cell_vertex(i);
+    return Quadrature<dim>(points);


Isn't this what QTrapez already does?

I'm wondering about this because the code you have here creates the same quadrature object every time. This might not matter because we don't create MappingFEField very often, but it's still annoying. It would be nice if you could just write something like

template <int dim> const Quadrature<dim> & // note: reference! get_vertex_quadrature () { static const QTrapez<dim> q; return q; }

which would create the quadrature formula only once.

QTrapez does not assure that the order of the quadrature points is the same as the order of the vertices, or does it? I'll make the code return a reference to a static object.

bangerth · 2018-04-25T03:51:15Z

source/fe/mapping_fe_field.cc

 inline
 const double &
-MappingFEField<dim,spacedim,DoFHandlerType,VectorType>::InternalData::shape
+MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData::shape


ha, how funny! nobody ever noticed this!

bangerth · 2018-04-25T03:51:40Z

source/fe/mapping_fe_field.cc

+  // we transform our tria iterator into a dof iterator so we can access
+  // data not associated with triangulations
+  typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell,
+                                                            euler_dof_handler);


make it const

luca-heltai · 2018-04-25T08:59:09Z

@masterleinad , this problem is related to your recent changes on the fe collection inside dh:

161/3560 Test #163: mappings/mapping_fe_field_02.debug ...................................***Failed   17.64 sec
Test mappings/mapping_fe_field_02.debug: RUN
===============================   OUTPUT BEGIN  ===============================
[ 50%] Generating mapping_fe_field_02.debug/interrupt_guard.cc
Scanning dependencies of target mapping_fe_field_02.debug
[ 50%] Building CXX object CMakeFiles/mapping_fe_field_02.debug.dir/mapping_fe_field_02.cc.o
Building CXX object CMakeFiles/mapping_fe_field_02.debug.dir/mapping_fe_field_02.debug/interrupt_guard.cc.o
Linking CXX executable mapping_fe_field_02.debug/mapping_fe_field_02.debug
[100%] Built target mapping_fe_field_02.debug
Scanning dependencies of target mapping_fe_field_02.debug.diff
[100%] Generating mapping_fe_field_02.debug/output
mappings/mapping_fe_field_02.debug: BUILD successful.
mappings/mapping_fe_field_02.debug: RUN failed. ------ Return code 134
mappings/mapping_fe_field_02.debug: RUN failed. ------ Result: /home/bob/build-gcc/tests/mappings/mapping_fe_field_02.debug/failing_output
mappings/mapping_fe_field_02.debug: RUN failed. ------ Partial output:
JobId e68cde002913 Tue Apr 24 20:26:11 2018
DEAL::dim, spacedim: 1, 1
DEAL::cells: 1, dofs: 2

mappings/mapping_fe_field_02.debug: RUN failed. ------ Additional output on stdout/stderr:


--------------------------------------------------------
An error occurred in line <104> of file </home/bob/source/source/base/subscriptor.cc> in function
    void dealii::Subscriptor::check_no_subscribers() const
The violated condition was: 
    counter == 0
Additional information: 
    (none)

Do you still make sure that the fe inside the dof handler is not destroyed and regenerated if it has not changed?

The following simple test breaks again MappingFEField and MappingQEulerian (even if I did not touch MappingQEulerian):

  Triangulation<dim,spacedim> tria;
  GridGenerator::hyper_cube (tria);
  FE_Q<dim,spacedim> fe(1);

  DoFHandler<dim,spacedim> dh(tria);
  dh.distribute_dofs(fe);

  ... build one of MappingFEField, or MappingQEulerian here

  // This is the root of the cause (which is now internally in both MappingFEField and MappingQEulerian)
  // They both have an FEValues object, that stores a SmartPointer to the fe.
  // SmartPointer<const FiniteElement<dim,spacedim>> fe_p(&dh.get_fe());

  tria.refine_global (1);
  dh.distribute_dofs(fe);

the second distribute_dofs fails, because there is an FEValues object inside both MappingQEulerian and MappingFEField that stores a SmartPointer to fe.

Apparently the failing part is here:

 if (fe_collection == nullptr || &((*fe_collection)[0]) != &ff)
     fe_collection = std_cxx14::make_unique<hp::FECollection<dim, spacedim>>(ff);  // Triggers.

Somehow the fe_collection is being recreated even if it should not, in principle, since the finite element is not changed. Is this maybe related to your recent changes @masterleinad (#6305 ?)

luca-heltai · 2018-04-25T18:43:45Z

PR #6330 should fix the failing test in this PR too. I'll rebase once that is merged.

tamiko · 2018-04-27T13:24:57Z

Blocks #6156

tamiko · 2018-04-27T14:38:41Z

@luca-heltai Please rebase and fix the indenting. :-)

This test fails because of dealii#6294.

luca-heltai added WIP Do not merge ☠️ labels Apr 23, 2018

luca-heltai added this to the Release 9.0 milestone Apr 23, 2018

bangerth reviewed Apr 23, 2018

View reviewed changes

davydden reviewed Apr 23, 2018

View reviewed changes

luca-heltai force-pushed the failing-test branch from 85b9604 to 5f69c2a Compare April 24, 2018 09:08

luca-heltai added Mapping and removed Do not merge ☠️ WIP labels Apr 24, 2018

davydden reviewed Apr 24, 2018

View reviewed changes

davydden approved these changes Apr 24, 2018

View reviewed changes

bangerth reviewed Apr 25, 2018

View reviewed changes

luca-heltai mentioned this pull request Apr 25, 2018

Failing test for multiple distribute dofs. #6326

Closed

luca-heltai force-pushed the failing-test branch 3 times, most recently from 3799004 to 3d30e10 Compare April 25, 2018 18:43

luca-heltai mentioned this pull request Apr 25, 2018

Bug in DoFHandler and hp::DoFHandler #6325

Closed

luca-heltai added 3 commits April 27, 2018 20:53

Test for min/max diameter and FEFieldFunction.

0b20a8d

This test fails because of dealii#6294.

Added implementation of MappingFEField::get_vertices

c35f044

Changelog.

08122ba

luca-heltai force-pushed the failing-test branch from 3d30e10 to 08122ba Compare April 27, 2018 18:55

masterleinad merged commit 9a48660 into dealii:master Apr 30, 2018

luca-heltai deleted the failing-test branch April 30, 2018 11:41

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix bug in MappingFEField #6296

Fix bug in MappingFEField #6296

luca-heltai commented Apr 23, 2018 •

edited

bangerth Apr 23, 2018

luca-heltai Apr 23, 2018

luca-heltai Apr 23, 2018

davydden Apr 23, 2018

davydden Apr 23, 2018 •

edited

davydden Apr 23, 2018

davydden Apr 23, 2018

davydden Apr 24, 2018

luca-heltai Apr 24, 2018

davydden Apr 24, 2018

davydden Apr 24, 2018

jppelteret Apr 24, 2018

davydden Apr 24, 2018

jppelteret Apr 24, 2018

davydden Apr 24, 2018

davydden Apr 24, 2018

luca-heltai Apr 24, 2018

luca-heltai Apr 24, 2018

davydden Apr 24, 2018

davydden Apr 24, 2018

davydden commented Apr 24, 2018

davydden commented Apr 24, 2018

bangerth commented Apr 25, 2018

bangerth Apr 25, 2018

luca-heltai Apr 25, 2018

bangerth Apr 25, 2018

bangerth Apr 25, 2018

luca-heltai commented Apr 25, 2018 •

edited

luca-heltai commented Apr 25, 2018 •

edited

tamiko commented Apr 27, 2018

tamiko commented Apr 27, 2018

Fix bug in MappingFEField #6296

Fix bug in MappingFEField #6296

Conversation

luca-heltai commented Apr 23, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden Apr 23, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden commented Apr 24, 2018

davydden commented Apr 24, 2018

bangerth commented Apr 25, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luca-heltai commented Apr 25, 2018 • edited

luca-heltai commented Apr 25, 2018 • edited

tamiko commented Apr 27, 2018

tamiko commented Apr 27, 2018

luca-heltai commented Apr 23, 2018 •

edited

davydden Apr 23, 2018 •

edited

luca-heltai commented Apr 25, 2018 •

edited

luca-heltai commented Apr 25, 2018 •

edited