Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
data_out.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
20 
21 #include <deal.II/fe/fe.h>
22 #include <deal.II/fe/fe_dgq.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping_q1.h>
25 
26 #include <deal.II/grid/tria.h>
28 
29 #include <deal.II/hp/dof_handler.h>
30 #include <deal.II/hp/fe_values.h>
31 
33 
34 #include <sstream>
35 
37 
38 
39 namespace internal
40 {
41  namespace DataOutImplementation
42  {
43  template <int dim, int spacedim>
45  const unsigned int n_datasets,
46  const unsigned int n_subdivisions,
47  const std::vector<unsigned int> &n_postprocessor_outputs,
48  const ::hp::MappingCollection<dim, spacedim> &mapping,
49  const std::vector<
50  std::shared_ptr<::hp::FECollection<dim, spacedim>>>
51  & finite_elements,
52  const UpdateFlags update_flags,
53  const std::vector<std::vector<unsigned int>> &cell_to_patch_index_map)
54  : ParallelDataBase<dim, spacedim>(n_datasets,
55  n_subdivisions,
56  n_postprocessor_outputs,
57  mapping,
58  finite_elements,
59  update_flags,
60  false)
61  , cell_to_patch_index_map(&cell_to_patch_index_map)
62  {}
63  } // namespace DataOutImplementation
64 } // namespace internal
65 
66 
67 
68 template <int dim, typename DoFHandlerType>
70 {
71  // For the moment, just call the existing virtual functions. This
72  // preserves backward compatibility. When these deprecated functions are
73  // removed, we'll just inline their functionality into the lambda
74  // functions created here.
75  set_cell_selection(
76  [this](const Triangulation<dim, spacedim> &) {
77  return this->first_locally_owned_cell();
78  },
79  [this](const Triangulation<dim, spacedim> &, const cell_iterator &cell) {
80  return this->next_locally_owned_cell(cell);
81  });
82 }
83 
84 
85 
86 template <int dim, typename DoFHandlerType>
87 void
89  const std::pair<cell_iterator, unsigned int> *cell_and_index,
90  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
91  DoFHandlerType::space_dimension>
92  & scratch_data,
93  const unsigned int n_subdivisions,
94  const CurvedCellRegion curved_cell_region)
95 {
96  // first create the output object that we will write into
97  ::DataOutBase::Patch<DoFHandlerType::dimension,
98  DoFHandlerType::space_dimension>
99  patch;
101  patch.reference_cell = cell_and_index->first->reference_cell();
102 
103  // initialize FEValues
104  scratch_data.reinit_all_fe_values(this->dof_data, cell_and_index->first);
105 
107  &fe_patch_values = scratch_data.get_present_fe_values(0);
108 
109  // set the vertices of the patch. if the mapping does not preserve locations
110  // (e.g. MappingQEulerian), we need to compute the offset of the vertex for
111  // the graphical output. Otherwise, we can just use the vertex info.
112  for (const unsigned int vertex : cell_and_index->first->vertex_indices())
113  if (scratch_data.mapping_collection[0].preserves_vertex_locations())
114  patch.vertices[vertex] = cell_and_index->first->vertex(vertex);
115  else
116  patch.vertices[vertex] =
117  fe_patch_values.get_mapping().transform_unit_to_real_cell(
118  cell_and_index->first,
120 
121  const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
122 
123  scratch_data.patch_values_scalar.solution_values.resize(n_q_points);
124  scratch_data.patch_values_scalar.solution_gradients.resize(n_q_points);
125  scratch_data.patch_values_scalar.solution_hessians.resize(n_q_points);
126  scratch_data.patch_values_system.solution_values.resize(n_q_points);
127  scratch_data.patch_values_system.solution_gradients.resize(n_q_points);
128  scratch_data.patch_values_system.solution_hessians.resize(n_q_points);
129 
130  for (unsigned int dataset = 0;
131  dataset < scratch_data.postprocessed_values.size();
132  ++dataset)
133  if (scratch_data.postprocessed_values[dataset].size() != 0)
134  scratch_data.postprocessed_values[dataset].resize(n_q_points);
135 
136  // First fill the geometric information for the patch: Where are the
137  // nodes in question located.
138  //
139  // Depending on the requested output of curved cells, if necessary
140  // append the quadrature points to the last rows of the patch.data
141  // member. This is the case if we want to produce curved cells at the
142  // boundary and this cell actually is at the boundary, or else if we
143  // want to produce curved cells everywhere
144  //
145  // note: a cell is *always* at the boundary if dim<spacedim
146  if (curved_cell_region == curved_inner_cells ||
147  (curved_cell_region == curved_boundary &&
148  (cell_and_index->first->at_boundary() ||
149  (DoFHandlerType::dimension != DoFHandlerType::space_dimension))) ||
150  (cell_and_index->first->reference_cell() !=
151  ReferenceCells::get_hypercube<dim>()))
152  {
153  Assert(patch.space_dim == DoFHandlerType::space_dimension,
154  ExcInternalError());
155 
156  // set the flag indicating that for this cell the points are
157  // explicitly given
158  patch.points_are_available = true;
159 
160  // then resize the patch.data member in order to have enough memory for
161  // the quadrature points as well, and copy the quadrature points there
162  const std::vector<Point<DoFHandlerType::space_dimension>> &q_points =
163  fe_patch_values.get_quadrature_points();
164 
165  patch.data.reinit(scratch_data.n_datasets +
166  DoFHandlerType::space_dimension,
167  n_q_points);
168  for (unsigned int i = 0; i < DoFHandlerType::space_dimension; ++i)
169  for (unsigned int q = 0; q < n_q_points; ++q)
170  patch.data(patch.data.size(0) - DoFHandlerType::space_dimension + i,
171  q) = q_points[q][i];
172  }
173  else
174  {
175  patch.data.reinit(scratch_data.n_datasets, n_q_points);
176  patch.points_are_available = false;
177  }
178 
179 
180  // Next fill the information we get from DoF data
181  if (scratch_data.n_datasets > 0)
182  {
183  // counter for data records
184  unsigned int offset = 0;
185 
186  // first fill dof_data
187  unsigned int dataset_number = 0;
188  for (const auto &dataset : this->dof_data)
189  {
190  const FEValuesBase<DoFHandlerType::dimension,
191  DoFHandlerType::space_dimension>
192  &this_fe_patch_values =
193  scratch_data.get_present_fe_values(dataset_number);
194  const unsigned int n_components =
195  this_fe_patch_values.get_fe().n_components();
196 
198  *postprocessor = dataset->postprocessor;
199 
200  if (postprocessor != nullptr)
201  {
202  // We have to postprocess the data, so determine, which fields
203  // have to be updated
204  const UpdateFlags update_flags =
205  postprocessor->get_needed_update_flags();
206 
207  if ((n_components == 1) &&
208  (dataset->is_complex_valued() == false))
209  {
210  // At each point there is only one component of value,
211  // gradient etc. Based on the 'if' statement above, we
212  // know that the solution is scalar and real-valued, so
213  // we do not need to worry about getting any imaginary
214  // components to the postprocessor, and we can safely
215  // call the function that evaluates a scalar field
216  if (update_flags & update_values)
217  dataset->get_function_values(
218  this_fe_patch_values,
220  real_part,
221  scratch_data.patch_values_scalar.solution_values);
222 
223  if (update_flags & update_gradients)
224  dataset->get_function_gradients(
225  this_fe_patch_values,
227  real_part,
228  scratch_data.patch_values_scalar.solution_gradients);
229 
230  if (update_flags & update_hessians)
231  dataset->get_function_hessians(
232  this_fe_patch_values,
234  real_part,
235  scratch_data.patch_values_scalar.solution_hessians);
236 
237  // Also fill some of the other fields postprocessors may
238  // want to access.
239  if (update_flags & update_quadrature_points)
240  scratch_data.patch_values_scalar.evaluation_points =
241  this_fe_patch_values.get_quadrature_points();
242 
243  const typename DoFHandlerType::cell_iterator dh_cell(
244  &cell_and_index->first->get_triangulation(),
245  cell_and_index->first->level(),
246  cell_and_index->first->index(),
247  dataset->dof_handler);
248  scratch_data.patch_values_scalar.template set_cell<dim>(
249  dh_cell);
250 
251  // Finally call the postprocessor's function that
252  // deals with scalar inputs.
253  postprocessor->evaluate_scalar_field(
254  scratch_data.patch_values_scalar,
255  scratch_data.postprocessed_values[dataset_number]);
256  }
257  else
258  {
259  // At each point we now have to evaluate a vector valued
260  // function and its derivatives. It may be that the solution
261  // is scalar and complex-valued, but we treat this as a vector
262  // field with two components.
263 
264  // At this point, we need to ask how we fill the fields that
265  // we want to pass on to the postprocessor. If the field in
266  // question is real-valued, we'll just extract the (only)
267  // real component from the solution fields
268  if (dataset->is_complex_valued() == false)
269  {
270  scratch_data.resize_system_vectors(n_components);
271 
272  if (update_flags & update_values)
273  dataset->get_function_values(
274  this_fe_patch_values,
276  real_part,
277  scratch_data.patch_values_system.solution_values);
278 
279  if (update_flags & update_gradients)
280  dataset->get_function_gradients(
281  this_fe_patch_values,
283  real_part,
284  scratch_data.patch_values_system.solution_gradients);
285 
286  if (update_flags & update_hessians)
287  dataset->get_function_hessians(
288  this_fe_patch_values,
290  real_part,
291  scratch_data.patch_values_system.solution_hessians);
292  }
293  else
294  {
295  // The solution is complex-valued. Let's cover the scalar
296  // case first (i.e., one scalar but complex-valued field,
297  // which we will have to split into its real and imaginar
298  // parts).
299  if (n_components == 1)
300  {
301  scratch_data.resize_system_vectors(2);
302 
303  // First get the real component of the scalar solution
304  // and copy the data into the
305  // scratch_data.patch_values_system output fields
306  if (update_flags & update_values)
307  {
308  dataset->get_function_values(
309  this_fe_patch_values,
312  scratch_data.patch_values_scalar
313  .solution_values);
314 
315  for (unsigned int i = 0;
316  i < scratch_data.patch_values_scalar
317  .solution_values.size();
318  ++i)
319  {
321  scratch_data.patch_values_system
322  .solution_values[i]
323  .size(),
324  2);
325  scratch_data.patch_values_system
326  .solution_values[i][0] =
327  scratch_data.patch_values_scalar
328  .solution_values[i];
329  }
330  }
331 
332  if (update_flags & update_gradients)
333  {
334  dataset->get_function_gradients(
335  this_fe_patch_values,
338  scratch_data.patch_values_scalar
339  .solution_gradients);
340 
341  for (unsigned int i = 0;
342  i < scratch_data.patch_values_scalar
343  .solution_gradients.size();
344  ++i)
345  {
347  scratch_data.patch_values_system
348  .solution_values[i]
349  .size(),
350  2);
351  scratch_data.patch_values_system
352  .solution_gradients[i][0] =
353  scratch_data.patch_values_scalar
354  .solution_gradients[i];
355  }
356  }
357 
358  if (update_flags & update_hessians)
359  {
360  dataset->get_function_hessians(
361  this_fe_patch_values,
364  scratch_data.patch_values_scalar
365  .solution_hessians);
366 
367  for (unsigned int i = 0;
368  i < scratch_data.patch_values_scalar
369  .solution_hessians.size();
370  ++i)
371  {
373  scratch_data.patch_values_system
374  .solution_hessians[i]
375  .size(),
376  2);
377  scratch_data.patch_values_system
378  .solution_hessians[i][0] =
379  scratch_data.patch_values_scalar
380  .solution_hessians[i];
381  }
382  }
383 
384  // Now we also have to get the imaginary
385  // component of the scalar solution
386  // and copy the data into the
387  // scratch_data.patch_values_system output fields
388  // that follow the real one
389  if (update_flags & update_values)
390  {
391  dataset->get_function_values(
392  this_fe_patch_values,
395  scratch_data.patch_values_scalar
396  .solution_values);
397 
398  for (unsigned int i = 0;
399  i < scratch_data.patch_values_scalar
400  .solution_values.size();
401  ++i)
402  {
404  scratch_data.patch_values_system
405  .solution_values[i]
406  .size(),
407  2);
408  scratch_data.patch_values_system
409  .solution_values[i][1] =
410  scratch_data.patch_values_scalar
411  .solution_values[i];
412  }
413  }
414 
415  if (update_flags & update_gradients)
416  {
417  dataset->get_function_gradients(
418  this_fe_patch_values,
421  scratch_data.patch_values_scalar
422  .solution_gradients);
423 
424  for (unsigned int i = 0;
425  i < scratch_data.patch_values_scalar
426  .solution_gradients.size();
427  ++i)
428  {
430  scratch_data.patch_values_system
431  .solution_values[i]
432  .size(),
433  2);
434  scratch_data.patch_values_system
435  .solution_gradients[i][1] =
436  scratch_data.patch_values_scalar
437  .solution_gradients[i];
438  }
439  }
440 
441  if (update_flags & update_hessians)
442  {
443  dataset->get_function_hessians(
444  this_fe_patch_values,
447  scratch_data.patch_values_scalar
448  .solution_hessians);
449 
450  for (unsigned int i = 0;
451  i < scratch_data.patch_values_scalar
452  .solution_hessians.size();
453  ++i)
454  {
456  scratch_data.patch_values_system
457  .solution_hessians[i]
458  .size(),
459  2);
460  scratch_data.patch_values_system
461  .solution_hessians[i][1] =
462  scratch_data.patch_values_scalar
463  .solution_hessians[i];
464  }
465  }
466  }
467  else
468  {
469  scratch_data.resize_system_vectors(2 * n_components);
470 
471  // This is the vector-valued, complex-valued case. In
472  // essence, we just need to do the same as above,
473  // i.e., call the functions in dataset
474  // to retrieve first the real and then the imaginary
475  // part of the solution, then copy them to the
476  // scratch_data.patch_values_system. The difference to
477  // the scalar case is that there, we could (ab)use the
478  // scratch_data.patch_values_scalar members for first
479  // the real part and then the imaginary part, copying
480  // them into the scratch_data.patch_values_system
481  // variable one after the other. We can't do this
482  // here because the solution is vector-valued, and
483  // so using the single
484  // scratch_data.patch_values_system object doesn't
485  // work.
486  //
487  // Rather, we need to come up with a temporary object
488  // for this (one for each the values, the gradients,
489  // and the hessians).
490  //
491  // Compared to the previous code path, we here
492  // first get the real and imaginary parts of the
493  // values, then the real and imaginary parts of the
494  // gradients, etc. This allows us to scope the
495  // temporary objects better
496  if (update_flags & update_values)
497  {
498  std::vector<Vector<double>> tmp(
499  scratch_data.patch_values_system.solution_values
500  .size(),
501  Vector<double>(n_components));
502 
503  // First get the real part into the tmp object
504  dataset->get_function_values(
505  this_fe_patch_values,
508  tmp);
509 
510  // Then copy these values into the first
511  // n_components slots of the output object.
512  for (unsigned int i = 0;
513  i < scratch_data.patch_values_system
514  .solution_values.size();
515  ++i)
516  {
518  scratch_data.patch_values_system
519  .solution_values[i]
520  .size(),
521  2 * n_components);
522  for (unsigned int j = 0; j < n_components;
523  ++j)
524  scratch_data.patch_values_system
525  .solution_values[i][j] = tmp[i][j];
526  }
527 
528  // Then do the same with the imaginary part,
529  // copying past the end of the previous set of
530  // values.
531  dataset->get_function_values(
532  this_fe_patch_values,
535  tmp);
536 
537  for (unsigned int i = 0;
538  i < scratch_data.patch_values_system
539  .solution_values.size();
540  ++i)
541  {
542  for (unsigned int j = 0; j < n_components;
543  ++j)
544  scratch_data.patch_values_system
545  .solution_values[i][j + n_components] =
546  tmp[i][j];
547  }
548  }
549 
550  // Now do the exact same thing for the gradients
551  if (update_flags & update_gradients)
552  {
553  std::vector<std::vector<
555  tmp(
556  scratch_data.patch_values_system
557  .solution_gradients.size(),
558  std::vector<
560  n_components));
561 
562  // First the real part
563  dataset->get_function_gradients(
564  this_fe_patch_values,
567  tmp);
568 
569  for (unsigned int i = 0;
570  i < scratch_data.patch_values_system
571  .solution_gradients.size();
572  ++i)
573  {
575  scratch_data.patch_values_system
576  .solution_gradients[i]
577  .size(),
578  2 * n_components);
579  for (unsigned int j = 0; j < n_components;
580  ++j)
581  scratch_data.patch_values_system
582  .solution_gradients[i][j] = tmp[i][j];
583  }
584 
585  // Then the imaginary part
586  dataset->get_function_gradients(
587  this_fe_patch_values,
590  tmp);
591 
592  for (unsigned int i = 0;
593  i < scratch_data.patch_values_system
594  .solution_gradients.size();
595  ++i)
596  {
597  for (unsigned int j = 0; j < n_components;
598  ++j)
599  scratch_data.patch_values_system
600  .solution_gradients[i][j + n_components] =
601  tmp[i][j];
602  }
603  }
604 
605  // And finally the Hessians. Same scheme as above.
606  if (update_flags & update_hessians)
607  {
608  std::vector<std::vector<
610  tmp(
611  scratch_data.patch_values_system
612  .solution_gradients.size(),
613  std::vector<
615  n_components));
616 
617  // First the real part
618  dataset->get_function_hessians(
619  this_fe_patch_values,
622  tmp);
623 
624  for (unsigned int i = 0;
625  i < scratch_data.patch_values_system
626  .solution_hessians.size();
627  ++i)
628  {
630  scratch_data.patch_values_system
631  .solution_hessians[i]
632  .size(),
633  2 * n_components);
634  for (unsigned int j = 0; j < n_components;
635  ++j)
636  scratch_data.patch_values_system
637  .solution_hessians[i][j] = tmp[i][j];
638  }
639 
640  // Then the imaginary part
641  dataset->get_function_hessians(
642  this_fe_patch_values,
645  tmp);
646 
647  for (unsigned int i = 0;
648  i < scratch_data.patch_values_system
649  .solution_hessians.size();
650  ++i)
651  {
652  for (unsigned int j = 0; j < n_components;
653  ++j)
654  scratch_data.patch_values_system
655  .solution_hessians[i][j + n_components] =
656  tmp[i][j];
657  }
658  }
659  }
660  }
661 
662  // Now set other fields we may need
663  if (update_flags & update_quadrature_points)
664  scratch_data.patch_values_system.evaluation_points =
665  this_fe_patch_values.get_quadrature_points();
666 
667  const typename DoFHandlerType::cell_iterator dh_cell(
668  &cell_and_index->first->get_triangulation(),
669  cell_and_index->first->level(),
670  cell_and_index->first->index(),
671  dataset->dof_handler);
672  scratch_data.patch_values_system.template set_cell<dim>(
673  dh_cell);
674 
675  // Whether the solution was complex-scalar or
676  // complex-vector-valued doesn't matter -- we took it apart
677  // into several fields and so we have to call the
678  // evaluate_vector_field() function.
679  postprocessor->evaluate_vector_field(
680  scratch_data.patch_values_system,
681  scratch_data.postprocessed_values[dataset_number]);
682  }
683 
684  // Now we need to copy the result of the postprocessor to
685  // the Patch object where it can then be further processed
686  // by the functions in DataOutBase
687  for (unsigned int q = 0; q < n_q_points; ++q)
688  for (unsigned int component = 0;
689  component < dataset->n_output_variables;
690  ++component)
691  patch.data(offset + component, q) =
692  scratch_data.postprocessed_values[dataset_number][q](
693  component);
694 
695  // Move the counter for the output location forward as
696  // appropriate
697  offset += dataset->n_output_variables;
698  }
699  else
700  {
701  // use the given data vector directly, without a postprocessor.
702  // again, we treat single component functions separately for
703  // efficiency reasons.
704  if (n_components == 1)
705  {
706  Assert(dataset->n_output_variables == 1, ExcInternalError());
707 
708  // First output the real part of the solution vector
709  dataset->get_function_values(
710  this_fe_patch_values,
712  real_part,
713  scratch_data.patch_values_scalar.solution_values);
714  for (unsigned int q = 0; q < n_q_points; ++q)
715  patch.data(offset, q) =
716  scratch_data.patch_values_scalar.solution_values[q];
717  offset += 1;
718 
719  // And if there is one, also output the imaginary part. Note
720  // that the problem is scalar-valued, so we can freely add the
721  // imaginary part after the real part without having to worry
722  // that we are interleaving the real components of a vector
723  // with the imaginary components of the same vector.
724  if (dataset->is_complex_valued() == true)
725  {
726  dataset->get_function_values(
727  this_fe_patch_values,
730  scratch_data.patch_values_scalar.solution_values);
731  for (unsigned int q = 0; q < n_q_points; ++q)
732  patch.data(offset, q) =
733  scratch_data.patch_values_scalar.solution_values[q];
734  offset += 1;
735  }
736  }
737  else
738  {
739  scratch_data.resize_system_vectors(n_components);
740 
741  // So we have a multi-component DoFHandler here. That's more
742  // complicated. If the vector is real-valued, then we can just
743  // get everything at all quadrature points and copy them into
744  // the output array. In fact, we don't have to worry at all
745  // about the interpretation of the components.
746  if (dataset->is_complex_valued() == false)
747  {
748  dataset->get_function_values(
749  this_fe_patch_values,
751  real_part,
752  scratch_data.patch_values_system.solution_values);
753  for (unsigned int component = 0; component < n_components;
754  ++component)
755  for (unsigned int q = 0; q < n_q_points; ++q)
756  patch.data(offset + component, q) =
757  scratch_data.patch_values_system.solution_values[q](
758  component);
759 
760  // Increment the counter for the actual data record.
761  offset += dataset->n_output_variables;
762  }
763  else
764  // The situation is more complicated if the input vector is
765  // complex-valued. The easiest approach would have been to
766  // just have all real and then all imaginary components.
767  // This would have been conceptually easy, but it has the
768  // annoying downside that if you have a vector-valued
769  // problem (say, [u v]) then the output order would have
770  // been [u_re, v_re, u_im, v_im]. That's tolerable, but not
771  // quite so nice because one typically thinks of real and
772  // imaginary parts as belonging together. We would really
773  // like the output order to be [u_re, u_im, v_re, v_im].
774  // That, too, would have been easy to implement because one
775  // just has to interleave real and imaginary parts.
776  //
777  // But that's also not what we want. That's because if one
778  // were, for example, to solve a complex-valued Stokes
779  // problem (e.g., computing eigenfunctions of the Stokes
780  // operator), then one has solution components
781  // [[u v] p] and the proper output order is
782  // [[u_re v_re] [u_im v_im] p_re p_im].
783  // In other words, the order in which we want to output
784  // data depends on the *interpretation* of components.
785  //
786  // Doing this requires a bit more code, and also needs to
787  // be in sync with what we do in
788  // DataOut_DoFData::get_dataset_names() and
789  // DataOut_DoFData::get_nonscalar_data_ranges().
790  {
791  // Given this description, first get the real parts of
792  // all components:
793  dataset->get_function_values(
794  this_fe_patch_values,
796  real_part,
797  scratch_data.patch_values_system.solution_values);
798 
799  // Then we need to distribute them to the correct
800  // location. This requires knowledge of the interpretation
801  // of components as discussed above.
802  {
803  Assert(dataset->data_component_interpretation.size() ==
804  n_components,
805  ExcInternalError());
806 
807  unsigned int destination = offset;
808  for (unsigned int component = 0;
809  component < n_components;
810  /* component is updated below */)
811  {
812  switch (
813  dataset->data_component_interpretation[component])
814  {
817  {
818  // OK, a scalar component. Put all of the
819  // values into the current row
820  // ('destination'); then move 'component'
821  // forward by one (so we treat the next
822  // component) and 'destination' forward by
823  // two (because we're going to put the
824  // imaginary part of the current component
825  // into the next slot).
826  for (unsigned int q = 0; q < n_q_points;
827  ++q)
828  patch.data(destination, q) =
829  scratch_data.patch_values_system
830  .solution_values[q](component);
831 
832  ++component;
833  destination += 2;
834 
835  break;
836  }
837 
840  {
841  // A vector component. Put the
842  // DoFHandlerType::space_dimension
843  // components into the next set of
844  // contiguous rows
845  // ('destination+c'); then move 'component'
846  // forward by spacedim (so we get to the
847  // next component after the current vector)
848  // and 'destination' forward by two*spacedim
849  // (because we're going to put the imaginary
850  // part of the vector into the subsequent
851  // spacedim slots).
852  const unsigned int size =
853  DoFHandlerType::space_dimension;
854  for (unsigned int c = 0; c < size; ++c)
855  for (unsigned int q = 0; q < n_q_points;
856  ++q)
857  patch.data(destination + c, q) =
858  scratch_data.patch_values_system
859  .solution_values[q](component + c);
860 
861  component += size;
862  destination += 2 * size;
863 
864  break;
865  }
866 
869  {
870  // Same approach as for vectors above.
871  const unsigned int size =
872  DoFHandlerType::space_dimension *
873  DoFHandlerType::space_dimension;
874  for (unsigned int c = 0; c < size; ++c)
875  for (unsigned int q = 0; q < n_q_points;
876  ++q)
877  patch.data(destination + c, q) =
878  scratch_data.patch_values_system
879  .solution_values[q](component + c);
880 
881  component += size;
882  destination += 2 * size;
883 
884  break;
885  }
886 
887  default:
888  Assert(false, ExcNotImplemented());
889  }
890  }
891  }
892 
893  // And now we need to do the same thing again for the
894  // imaginary parts, starting at the top of the list of
895  // components/destinations again.
896  dataset->get_function_values(
897  this_fe_patch_values,
900  scratch_data.patch_values_system.solution_values);
901  {
902  unsigned int destination = offset;
903  for (unsigned int component = 0;
904  component < n_components;
905  /* component is updated below */)
906  {
907  switch (
908  dataset->data_component_interpretation[component])
909  {
912  {
913  // OK, a scalar component. Put all of the
914  // values into the row past the current one
915  // ('destination+1') since 'destination' is
916  // occupied by the real part.
917  for (unsigned int q = 0; q < n_q_points;
918  ++q)
919  patch.data(destination + 1, q) =
920  scratch_data.patch_values_system
921  .solution_values[q](component);
922 
923  ++component;
924  destination += 2;
925 
926  break;
927  }
928 
931  {
932  // A vector component. Put the
933  // DoFHandlerType::space_dimension
934  // components into the set of contiguous
935  // rows that follow the real parts
936  // ('destination+spacedim+c').
937  const unsigned int size =
938  DoFHandlerType::space_dimension;
939  for (unsigned int c = 0; c < size; ++c)
940  for (unsigned int q = 0; q < n_q_points;
941  ++q)
942  patch.data(destination + size + c, q) =
943  scratch_data.patch_values_system
944  .solution_values[q](component + c);
945 
946  component += size;
947  destination += 2 * size;
948 
949  break;
950  }
951 
954  {
955  // Same as for vectors.
956  const unsigned int size =
957  DoFHandlerType::space_dimension *
958  DoFHandlerType::space_dimension;
959  for (unsigned int c = 0; c < size; ++c)
960  for (unsigned int q = 0; q < n_q_points;
961  ++q)
962  patch.data(destination + size + c, q) =
963  scratch_data.patch_values_system
964  .solution_values[q](component + c);
965 
966  component += size;
967  destination += 2 * size;
968 
969  break;
970  }
971 
972  default:
973  Assert(false, ExcNotImplemented());
974  }
975  }
976  }
977 
978  // Increment the counter for the actual data record. We
979  // need to move it forward a number of positions equal to
980  // the number of components of this data set, times two
981  // because we dealt with a complex-valued input vector
982  offset += dataset->n_output_variables * 2;
983  }
984  }
985  }
986 
987  // Also update the dataset_number index that we carry along with the
988  // for-loop over all data sets.
989  ++dataset_number;
990  }
991 
992  // Then do the cell data. At least, we don't have to worry about
993  // complex-valued vectors/tensors since cell data is always scalar.
994  if (this->cell_data.size() != 0)
995  {
996  Assert(!cell_and_index->first->has_children(), ExcNotImplemented());
997 
998  for (const auto &dataset : this->cell_data)
999  {
1000  // as above, first output the real part
1001  {
1002  const double value =
1003  dataset->get_cell_data_value(cell_and_index->second,
1006  for (unsigned int q = 0; q < n_q_points; ++q)
1007  patch.data(offset, q) = value;
1008  }
1009 
1010  // and if there is one, also output the imaginary part
1011  if (dataset->is_complex_valued() == true)
1012  {
1013  const double value = dataset->get_cell_data_value(
1014  cell_and_index->second,
1016  imaginary_part);
1017  for (unsigned int q = 0; q < n_q_points; ++q)
1018  patch.data(offset + 1, q) = value;
1019  }
1020 
1021  offset += (dataset->is_complex_valued() ? 2 : 1);
1022  }
1023  }
1024  }
1025 
1026 
1027  for (const unsigned int f : cell_and_index->first->face_indices())
1028  {
1029  // let's look up whether the neighbor behind that face is noted in the
1030  // table of cells which we treat. this can only happen if the neighbor
1031  // exists, and is on the same level as this cell, but it may also happen
1032  // that the neighbor is not a member of the range of cells over which we
1033  // loop, in which case the respective entry in the
1034  // cell_to_patch_index_map will have the value no_neighbor. (note that
1035  // since we allocated only as much space in this array as the maximum
1036  // index of the cells we loop over, not every neighbor may have its
1037  // space in it, so we have to assume that it is extended by values
1038  // no_neighbor)
1039  if (cell_and_index->first->at_boundary(f) ||
1040  (cell_and_index->first->neighbor(f)->level() !=
1041  cell_and_index->first->level()))
1042  {
1043  patch.neighbors[f] = numbers::invalid_unsigned_int;
1044  continue;
1045  }
1046 
1047  const cell_iterator neighbor = cell_and_index->first->neighbor(f);
1048  Assert(static_cast<unsigned int>(neighbor->level()) <
1049  scratch_data.cell_to_patch_index_map->size(),
1050  ExcInternalError());
1051  if ((static_cast<unsigned int>(neighbor->index()) >=
1052  (*scratch_data.cell_to_patch_index_map)[neighbor->level()].size()) ||
1053  ((*scratch_data.cell_to_patch_index_map)[neighbor->level()]
1054  [neighbor->index()] ==
1056  {
1057  patch.neighbors[f] = numbers::invalid_unsigned_int;
1058  continue;
1059  }
1060 
1061  // now, there is a neighbor, so get its patch number and set it for the
1062  // neighbor index
1063  patch.neighbors[f] =
1064  (*scratch_data
1065  .cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
1066  }
1067 
1068  const unsigned int patch_idx =
1069  (*scratch_data.cell_to_patch_index_map)[cell_and_index->first->level()]
1070  [cell_and_index->first->index()];
1071  // did we mess up the indices?
1072  Assert(patch_idx < this->patches.size(), ExcInternalError());
1073  patch.patch_index = patch_idx;
1074 
1075  // Put the patch into the patches vector. instead of copying the data,
1076  // simply swap the contents to avoid the penalty of writing into another
1077  // processor's memory
1078  this->patches[patch_idx].swap(patch);
1079 }
1080 
1081 
1082 
1083 template <int dim, typename DoFHandlerType>
1084 void
1086 {
1087  AssertDimension(this->triangulation->get_reference_cells().size(), 1);
1088 
1089  build_patches(
1090  this->triangulation->get_reference_cells()[0]
1091  .template get_default_linear_mapping<DoFHandlerType::dimension,
1092  DoFHandlerType::space_dimension>(),
1093  n_subdivisions,
1094  no_curved_cells);
1095 }
1096 
1097 
1098 
1099 template <int dim, typename DoFHandlerType>
1100 void
1103  & mapping,
1104  const unsigned int n_subdivisions_,
1105  const CurvedCellRegion curved_region)
1106 {
1107  hp::MappingCollection<DoFHandlerType::dimension,
1108  DoFHandlerType::space_dimension>
1109  mapping_collection(mapping);
1110 
1111  build_patches(mapping_collection, n_subdivisions_, curved_region);
1112 }
1113 
1114 
1115 
1116 template <int dim, typename DoFHandlerType>
1117 void
1119  const hp::MappingCollection<DoFHandlerType::dimension,
1120  DoFHandlerType::space_dimension> &mapping,
1121  const unsigned int n_subdivisions_,
1122  const CurvedCellRegion curved_region)
1123 {
1124  // Check consistency of redundant template parameter
1125  Assert(dim == DoFHandlerType::dimension,
1126  ExcDimensionMismatch(dim, DoFHandlerType::dimension));
1127 
1128  Assert(this->triangulation != nullptr,
1130 
1131  const unsigned int n_subdivisions =
1132  (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
1133  Assert(n_subdivisions >= 1,
1135  n_subdivisions));
1136 
1137  this->validate_dataset_names();
1138 
1139  // First count the cells we want to create patches of. Also fill the object
1140  // that maps the cell indices to the patch numbers, as this will be needed
1141  // for generation of neighborship information.
1142  // Note, there is a confusing mess of different indices here at play:
1143  // - patch_index: the index of a patch in all_cells
1144  // - cell->index: only unique on each level, used in cell_to_patch_index_map
1145  // - active_index: index for a cell when counting from begin_active() using
1146  // ++cell (identical to cell->active_cell_index())
1147  // - cell_index: unique index of a cell counted using
1148  // next_cell_function() starting from first_cell_function()
1149  //
1150  // It turns out that we create one patch for each selected cell, so
1151  // patch_index==cell_index.
1152  //
1153  // Now construct the map such that
1154  // cell_to_patch_index_map[cell->level][cell->index] = patch_index
1155  std::vector<std::vector<unsigned int>> cell_to_patch_index_map;
1156  cell_to_patch_index_map.resize(this->triangulation->n_levels());
1157  for (unsigned int l = 0; l < this->triangulation->n_levels(); ++l)
1158  {
1159  // max_index is the largest cell->index on level l
1160  unsigned int max_index = 0;
1161  for (cell_iterator cell = first_cell_function(*this->triangulation);
1162  cell != this->triangulation->end();
1163  cell = next_cell_function(*this->triangulation, cell))
1164  if (static_cast<unsigned int>(cell->level()) == l)
1165  max_index =
1166  std::max(max_index, static_cast<unsigned int>(cell->index()));
1167 
1168  cell_to_patch_index_map[l].resize(
1169  max_index + 1,
1171  DoFHandlerType::dimension,
1172  DoFHandlerType::space_dimension>::no_neighbor);
1173  }
1174 
1175  // will be all_cells[patch_index] = pair(cell, active_index)
1176  std::vector<std::pair<cell_iterator, unsigned int>> all_cells;
1177  {
1178  // important: we need to compute the active_index of the cell in the range
1179  // 0..n_active_cells() because this is where we need to look up cell
1180  // data from (cell data vectors do not have the length distance computed by
1181  // first_cell_function/next_cell_function because this might skip
1182  // some values (FilteredIterator).
1183  auto active_cell = this->triangulation->begin_active();
1184  unsigned int active_index = 0;
1185  cell_iterator cell = first_cell_function(*this->triangulation);
1186  for (; cell != this->triangulation->end();
1187  cell = next_cell_function(*this->triangulation, cell))
1188  {
1189  // move forward until active_cell points at the cell (cell) we are
1190  // looking at to compute the current active_index
1191  while (active_cell != this->triangulation->end() && cell->is_active() &&
1192  decltype(active_cell)(cell) != active_cell)
1193  {
1194  ++active_cell;
1195  ++active_index;
1196  }
1197 
1198  Assert(static_cast<unsigned int>(cell->level()) <
1199  cell_to_patch_index_map.size(),
1200  ExcInternalError());
1201  Assert(static_cast<unsigned int>(cell->index()) <
1202  cell_to_patch_index_map[cell->level()].size(),
1203  ExcInternalError());
1204  Assert(active_index < this->triangulation->n_active_cells(),
1205  ExcInternalError());
1206  cell_to_patch_index_map[cell->level()][cell->index()] =
1207  all_cells.size();
1208 
1209  all_cells.emplace_back(cell, active_index);
1210  }
1211  }
1212 
1213  this->patches.clear();
1214  this->patches.resize(all_cells.size());
1215 
1216  // Now create a default object for the WorkStream object to work with. The
1217  // first step is to count how many output data sets there will be. This is,
1218  // in principle, just the number of components of each data set, but we
1219  // need to allocate two entries per component if there are
1220  // complex-valued input data (unless we use a postprocessor on this
1221  // output -- all postprocessor outputs are real-valued)
1222  unsigned int n_datasets = 0;
1223  for (unsigned int i = 0; i < this->cell_data.size(); ++i)
1224  n_datasets += (this->cell_data[i]->is_complex_valued() &&
1225  (this->cell_data[i]->postprocessor == nullptr) ?
1226  2 :
1227  1);
1228  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
1229  n_datasets += (this->dof_data[i]->n_output_variables *
1230  (this->dof_data[i]->is_complex_valued() &&
1231  (this->dof_data[i]->postprocessor == nullptr) ?
1232  2 :
1233  1));
1234 
1235  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
1236  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
1237  if (this->dof_data[dataset]->postprocessor)
1238  n_postprocessor_outputs[dataset] =
1239  this->dof_data[dataset]->n_output_variables;
1240  else
1241  n_postprocessor_outputs[dataset] = 0;
1242 
1243  const CurvedCellRegion curved_cell_region =
1244  (n_subdivisions < 2 ? no_curved_cells : curved_region);
1245 
1247  if (curved_cell_region != no_curved_cells)
1248  update_flags |= update_quadrature_points;
1249 
1250  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
1251  if (this->dof_data[i]->postprocessor)
1252  update_flags |=
1253  this->dof_data[i]->postprocessor->get_needed_update_flags();
1254  // perhaps update_normal_vectors is present, which would only be useful on
1255  // faces, but we may not use it here.
1256  Assert(
1257  !(update_flags & update_normal_vectors),
1258  ExcMessage(
1259  "The update of normal vectors may not be requested for evaluation of "
1260  "data on cells via DataPostprocessor."));
1261 
1262  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
1263  DoFHandlerType::space_dimension>
1264  thread_data(n_datasets,
1265  n_subdivisions,
1266  n_postprocessor_outputs,
1267  mapping,
1268  this->get_fes(),
1269  update_flags,
1270  cell_to_patch_index_map);
1271 
1272  auto worker = [this, n_subdivisions, curved_cell_region](
1273  const std::pair<cell_iterator, unsigned int> *cell_and_index,
1275  DoFHandlerType::dimension,
1276  DoFHandlerType::space_dimension> &scratch_data,
1277  // this function doesn't actually need a copy data object --
1278  // it just writes everything right into the output array
1279  int) {
1280  this->build_one_patch(cell_and_index,
1281  scratch_data,
1282  n_subdivisions,
1283  curved_cell_region);
1284  };
1285 
1286  // now build the patches in parallel
1287  if (all_cells.size() > 0)
1288  WorkStream::run(all_cells.data(),
1289  all_cells.data() + all_cells.size(),
1290  worker,
1291  // no copy-local-to-global function needed here
1292  std::function<void(const int)>(),
1293  thread_data,
1294  /* dummy CopyData object = */ 0,
1295  // experimenting shows that we can make things run a bit
1296  // faster if we increase the number of cells we work on
1297  // per item (i.e., WorkStream's chunk_size argument,
1298  // about 10% improvement) and the items in flight at any
1299  // given time (another 5% on the testcase discussed in
1300  // @ref workstream_paper, on 32 cores) and if
1302  64);
1303 }
1304 
1305 
1306 
1307 template <int dim, typename DoFHandlerType>
1308 void
1310  const std::function<cell_iterator(const Triangulation<dim, spacedim> &)>
1311  & first_cell,
1312  const std::function<cell_iterator(const Triangulation<dim, spacedim> &,
1313  const cell_iterator &)> &next_cell)
1314 {
1315  first_cell_function = first_cell;
1316  next_cell_function = next_cell;
1317 }
1318 
1319 
1320 
1321 template <int dim, typename DoFHandlerType>
1322 void
1324  const FilteredIterator<cell_iterator> &filtered_iterator)
1325 {
1326  const auto first_cell =
1327  [filtered_iterator](const Triangulation<dim, spacedim> &triangulation) {
1328  // Create a copy of the filtered iterator so that we can
1329  // call a non-const function -- though we are really only
1330  // interested in the return value of that function, not the
1331  // state of the object
1332  FilteredIterator<cell_iterator> x = filtered_iterator;
1333  return x.set_to_next_positive(triangulation.begin());
1334  };
1335 
1336 
1337  const auto next_cell =
1338  [filtered_iterator](const Triangulation<dim, spacedim> &,
1339  const cell_iterator &cell) {
1340  // Create a copy of the filtered iterator so that we can
1341  // call a non-const function -- though we are really only
1342  // interested in the return value of that function, not the
1343  // state of the object
1344  FilteredIterator<cell_iterator> x = filtered_iterator;
1345 
1346  // Set the iterator to 'cell'. Since 'cell' must satisfy the
1347  // predicate (that's how it was created), set_to_next_positive
1348  // simply sets the iterator to 'cell'.
1349  x.set_to_next_positive(cell);
1350 
1351  // Advance by one:
1352  ++x;
1353 
1354  return x;
1355  };
1356 
1357  set_cell_selection(first_cell, next_cell);
1358 }
1359 
1360 
1361 
1362 template <int dim, typename DoFHandlerType>
1363 const std::pair<typename DataOut<dim, DoFHandlerType>::FirstCellFunctionType,
1366 {
1367  return std::make_pair(first_cell_function, next_cell_function);
1368 }
1369 
1370 
1371 
1372 template <int dim, typename DoFHandlerType>
1375 {
1376  return this->triangulation->begin_active();
1377 }
1378 
1379 
1380 
1381 template <int dim, typename DoFHandlerType>
1384  const typename DataOut<dim, DoFHandlerType>::cell_iterator &cell)
1385 {
1386  // convert the iterator to an active_iterator and advance this to the next
1387  // active cell
1388  typename Triangulation<DoFHandlerType::dimension,
1389  DoFHandlerType::space_dimension>::active_cell_iterator
1390  active_cell = cell;
1391  ++active_cell;
1392  return active_cell;
1393 }
1394 
1395 
1396 
1397 template <int dim, typename DoFHandlerType>
1400 {
1401  typename DataOut<dim, DoFHandlerType>::cell_iterator cell = first_cell();
1402 
1403  // skip cells if the current one has no children (is active) and is a ghost
1404  // or artificial cell
1405  while ((cell != this->triangulation->end()) && cell->is_active() &&
1406  !cell->is_locally_owned())
1407  cell = next_cell(cell);
1408 
1409  return cell;
1410 }
1411 
1412 
1413 
1414 template <int dim, typename DoFHandlerType>
1417  const typename DataOut<dim, DoFHandlerType>::cell_iterator &old_cell)
1418 {
1420  next_cell(old_cell);
1421  while ((cell != this->triangulation->end()) && cell->is_active() &&
1422  !cell->is_locally_owned())
1423  cell = next_cell(cell);
1424  return cell;
1425 }
1426 
1427 
1428 // explicit instantiations
1429 #include "data_out.inst"
1430 
Shape function values.
typename std::function< cell_iterator(const Triangulation< dim, spacedim > &, const cell_iterator &)> NextCellFunctionType
Definition: data_out.h:180
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static ::ExceptionBase & ExcNoTriangulationSelected()
ParallelData(const unsigned int n_datasets, const unsigned int n_subdivisions, const std::vector< unsigned int > &n_postprocessor_outputs, const ::hp::MappingCollection< dim, spacedim > &mapping, const std::vector< std::shared_ptr<::hp::FECollection< dim, spacedim >>> &finite_elements, const UpdateFlags update_flags, const std::vector< std::vector< unsigned int >> &cell_to_patch_index_map)
Definition: data_out.cc:44
void build_one_patch(const std::pair< cell_iterator, unsigned int > *cell_and_index, internal::DataOutImplementation::ParallelData< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &scratch_data, const unsigned int n_subdivisions, const CurvedCellRegion curved_cell_region)
Definition: data_out.cc:88
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
DataOut()
Definition: data_out.cc:69
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:240
const std::vector< std::vector< unsigned int > > * cell_to_patch_index_map
Definition: data_out.h:56
const Mapping< dim, spacedim > & get_mapping() const
virtual cell_iterator next_cell(const cell_iterator &cell)
Definition: data_out.cc:1383
typename DataOut_DoFData< DoFHandler< dim >, DoFHandler< dim > ::dimension, DoFHandler< dim > ::space_dimension >::cell_iterator cell_iterator
Definition: data_out.h:165
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Transformed quadrature points.
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
const std::pair< FirstCellFunctionType, NextCellFunctionType > get_cell_selection() const
Definition: data_out.cc:1365
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1085
const std::vector< Point< spacedim > > & get_quadrature_points() const
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual cell_iterator first_locally_owned_cell()
Definition: data_out.cc:1399
#define Assert(cond, exc)
Definition: exceptions.h:1465
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:303
virtual cell_iterator first_cell()
Definition: data_out.cc:1374
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
unsigned int n_subdivisions
Second derivatives of shape functions.
const ::hp::MappingCollection< dim, spacedim > mapping_collection
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
void set_cell_selection(const std::function< cell_iterator(const Triangulation< dim, spacedim > &)> &first_cell, const std::function< cell_iterator(const Triangulation< dim, spacedim > &, const cell_iterator &)> &next_cell)
Definition: data_out.cc:1309
FilteredIterator & set_to_next_positive(const BaseIterator &bi)
Definition: tensor.h:462
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
Shape function gradients.
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Normal vectors.
static ::ExceptionBase & ExcNotImplemented()
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1337
static unsigned int n_threads()
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
virtual cell_iterator next_locally_owned_cell(const cell_iterator &cell)
Definition: data_out.cc:1416
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
virtual UpdateFlags get_needed_update_flags() const =0