error_estimator.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifdef OOMPH_HAS_MPI
27 #include "mpi.h"
28 #endif
29 
30 
32 #include "error_estimator.h"
33 #include "shape.h"
34 #include "Telements.h"
35 
36 namespace oomph
37 {
38  //====================================================================
39  /// Recovery shape functions as functions of the global, Eulerian
40  /// coordinate x of dimension dim.
41  /// The recovery shape functions are complete polynomials of
42  /// the order specified by Recovery_order.
43  //====================================================================
45  const unsigned& dim,
46  Vector<double>& psi_r)
47  {
48  std::ostringstream error_stream;
49 
50  /// Which spatial dimension are we dealing with?
51  switch (dim)
52  {
53  case 1:
54 
55  // 1D:
56  //====
57 
58  /// Find order of recovery shape functions
59  switch (recovery_order())
60  {
61  case 1:
62 
63  // Complete linear polynomial in 1D:
64  psi_r[0] = 1.0;
65  psi_r[1] = x[0];
66  break;
67 
68  case 2:
69 
70  // Complete quadratic polynomial in 1D:
71  psi_r[0] = 1.0;
72  psi_r[1] = x[0];
73  psi_r[2] = x[0] * x[0];
74  break;
75 
76  case 3:
77 
78  // Complete cubic polynomial in 1D:
79  psi_r[0] = 1.0;
80  psi_r[1] = x[0];
81  psi_r[2] = x[0] * x[0];
82  psi_r[3] = x[0] * x[0] * x[0];
83  break;
84 
85  default:
86 
87  error_stream << "Recovery shape functions for recovery order "
88  << recovery_order()
89  << " haven't yet been implemented for 1D" << std::endl;
90 
91  throw OomphLibError(error_stream.str(),
92  OOMPH_CURRENT_FUNCTION,
93  OOMPH_EXCEPTION_LOCATION);
94  }
95  break;
96 
97  case 2:
98 
99  // 2D:
100  //====
101 
102  /// Find order of recovery shape functions
103  switch (recovery_order())
104  {
105  case 1:
106 
107  // Complete linear polynomial in 2D:
108  psi_r[0] = 1.0;
109  psi_r[1] = x[0];
110  psi_r[2] = x[1];
111  break;
112 
113  case 2:
114 
115  // Complete quadratic polynomial in 2D:
116  psi_r[0] = 1.0;
117  psi_r[1] = x[0];
118  psi_r[2] = x[1];
119  psi_r[3] = x[0] * x[0];
120  psi_r[4] = x[0] * x[1];
121  psi_r[5] = x[1] * x[1];
122  break;
123 
124  case 3:
125 
126  // Complete cubic polynomial in 2D:
127  psi_r[0] = 1.0;
128  psi_r[1] = x[0];
129  psi_r[2] = x[1];
130  psi_r[3] = x[0] * x[0];
131  psi_r[4] = x[0] * x[1];
132  psi_r[5] = x[1] * x[1];
133  psi_r[6] = x[0] * x[0] * x[0];
134  psi_r[7] = x[0] * x[0] * x[1];
135  psi_r[8] = x[0] * x[1] * x[1];
136  psi_r[9] = x[1] * x[1] * x[1];
137  break;
138 
139  default:
140 
141  error_stream << "Recovery shape functions for recovery order "
142  << recovery_order()
143  << " haven't yet been implemented for 2D" << std::endl;
144 
145  throw OomphLibError(error_stream.str(),
146  OOMPH_CURRENT_FUNCTION,
147  OOMPH_EXCEPTION_LOCATION);
148  }
149  break;
150 
151  case 3:
152 
153  // 3D:
154  //====
155  /// Find order of recovery shape functions
156  switch (recovery_order())
157  {
158  case 1:
159 
160  // Complete linear polynomial in 3D:
161  psi_r[0] = 1.0;
162  psi_r[1] = x[0];
163  psi_r[2] = x[1];
164  psi_r[3] = x[2];
165  break;
166 
167  case 2:
168 
169  // Complete quadratic polynomial in 3D:
170  psi_r[0] = 1.0;
171  psi_r[1] = x[0];
172  psi_r[2] = x[1];
173  psi_r[3] = x[2];
174  psi_r[4] = x[0] * x[0];
175  psi_r[5] = x[0] * x[1];
176  psi_r[6] = x[0] * x[2];
177  psi_r[7] = x[1] * x[1];
178  psi_r[8] = x[1] * x[2];
179  psi_r[9] = x[2] * x[2];
180  break;
181 
182  case 3:
183 
184  // Complete cubic polynomial in 3D:
185  psi_r[0] = 1.0;
186  psi_r[1] = x[0];
187  psi_r[2] = x[1];
188  psi_r[3] = x[2];
189  psi_r[4] = x[0] * x[0];
190  psi_r[5] = x[0] * x[1];
191  psi_r[6] = x[0] * x[2];
192  psi_r[7] = x[1] * x[1];
193  psi_r[8] = x[1] * x[2];
194  psi_r[9] = x[2] * x[2];
195  psi_r[10] = x[0] * x[0] * x[0];
196  psi_r[11] = x[0] * x[0] * x[1];
197  psi_r[12] = x[0] * x[0] * x[2];
198  psi_r[13] = x[1] * x[1] * x[1];
199  psi_r[14] = x[0] * x[1] * x[1];
200  psi_r[15] = x[2] * x[1] * x[1];
201  psi_r[16] = x[2] * x[2] * x[2];
202  psi_r[17] = x[2] * x[2] * x[0];
203  psi_r[18] = x[2] * x[2] * x[1];
204  psi_r[19] = x[0] * x[1] * x[2];
205 
206  break;
207 
208  default:
209 
210  error_stream << "Recovery shape functions for recovery order "
211  << recovery_order()
212  << " haven't yet been implemented for 3D" << std::endl;
213 
214  throw OomphLibError(error_stream.str(),
215  OOMPH_CURRENT_FUNCTION,
216  OOMPH_EXCEPTION_LOCATION);
217  }
218 
219 
220  break;
221 
222  default:
223 
224  // Any other dimension?
225  //=====================
226  error_stream << "No recovery shape functions for dim " << dim
227  << std::endl;
228  throw OomphLibError(
229  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
230  break;
231  }
232  }
233 
234  //====================================================================
235  /// Integation scheme associated with the recovery shape functions
236  /// must be of sufficiently high order to integrate the mass matrix
237  /// associated with the recovery shape functions. The argument
238  /// is the dimension of the elements.
239  /// The integration is performed locally over the elements, so the
240  /// integration scheme does depend on the geometry of the element.
241  /// The type of element is specified by the boolean which is
242  /// true if elements in the patch are QElements and false if they are
243  /// TElements (will need change if we ever have other element types)
244  //====================================================================
246  const bool& is_q_mesh)
247  {
248  std::ostringstream error_stream;
249 
250  /// Which spatial dimension are we dealing with?
251  switch (dim)
252  {
253  case 1:
254 
255  // 1D:
256  //====
257 
258  /// Find order of recovery shape functions
259  switch (recovery_order())
260  {
261  case 1:
262 
263  // Complete linear polynomial in 1D
264  //(quadratic terms in mass matrix)
265  if (is_q_mesh)
266  {
267  return (new Gauss<1, 2>);
268  }
269  else
270  {
271  return (new TGauss<1, 2>);
272  }
273  break;
274 
275  case 2:
276 
277  // Complete quadratic polynomial in 1D:
278  //(quartic terms in the mass marix)
279  if (is_q_mesh)
280  {
281  return (new Gauss<1, 3>);
282  }
283  else
284  {
285  return (new TGauss<1, 3>);
286  }
287  break;
288 
289  case 3:
290 
291  // Complete cubic polynomial in 1D:
292  // (order six terms in mass matrix)
293  if (is_q_mesh)
294  {
295  return (new Gauss<1, 4>);
296  }
297  else
298  {
299  return (new TGauss<1, 4>);
300  }
301  break;
302 
303  default:
304 
305  error_stream << "Recovery shape functions for recovery order "
306  << recovery_order()
307  << " haven't yet been implemented for 1D" << std::endl;
308 
309  throw OomphLibError(error_stream.str(),
310  OOMPH_CURRENT_FUNCTION,
311  OOMPH_EXCEPTION_LOCATION);
312  }
313  break;
314 
315  case 2:
316 
317  // 2D:
318  //====
319 
320  /// Find order of recovery shape functions
321  switch (recovery_order())
322  {
323  case 1:
324 
325  // Complete linear polynomial in 2D:
326  if (is_q_mesh)
327  {
328  return (new Gauss<2, 2>);
329  }
330  else
331  {
332  return (new TGauss<2, 2>);
333  }
334  break;
335 
336  case 2:
337 
338  // Complete quadratic polynomial in 2D:
339  if (is_q_mesh)
340  {
341  return (new Gauss<2, 3>);
342  }
343  else
344  {
345  return (new TGauss<2, 3>);
346  }
347  break;
348 
349  case 3:
350 
351  // Complete cubic polynomial in 2D:
352  if (is_q_mesh)
353  {
354  return (new Gauss<2, 4>);
355  }
356  else
357  {
358  return (new TGauss<2, 4>);
359  }
360  break;
361 
362  default:
363 
364  error_stream << "Recovery shape functions for recovery order "
365  << recovery_order()
366  << " haven't yet been implemented for 2D" << std::endl;
367 
368  throw OomphLibError(error_stream.str(),
369  OOMPH_CURRENT_FUNCTION,
370  OOMPH_EXCEPTION_LOCATION);
371  }
372  break;
373 
374  case 3:
375 
376  // 3D:
377  //====
378  /// Find order of recovery shape functions
379  switch (recovery_order())
380  {
381  case 1:
382 
383  // Complete linear polynomial in 3D:
384  if (is_q_mesh)
385  {
386  return (new Gauss<3, 2>);
387  }
388  else
389  {
390  return (new TGauss<3, 2>);
391  }
392  break;
393 
394  case 2:
395 
396  // Complete quadratic polynomial in 3D:
397  if (is_q_mesh)
398  {
399  return (new Gauss<3, 3>);
400  }
401  else
402  {
403  return (new TGauss<3, 3>);
404  }
405  break;
406 
407  case 3:
408 
409  // Complete cubic polynomial in 3D:
410  if (is_q_mesh)
411  {
412  return (new Gauss<3, 4>);
413  }
414  else
415  {
416  return (new TGauss<3, 5>);
417  } // TGauss<3,4> not implemented
418 
419  break;
420 
421  default:
422 
423  error_stream << "Recovery shape functions for recovery order "
424  << recovery_order()
425  << " haven't yet been implemented for 3D" << std::endl;
426 
427  throw OomphLibError(error_stream.str(),
428  OOMPH_CURRENT_FUNCTION,
429  OOMPH_EXCEPTION_LOCATION);
430  }
431 
432 
433  break;
434 
435  default:
436 
437  // Any other dimension?
438  //=====================
439  error_stream << "No recovery shape functions for dim " << dim
440  << std::endl;
441  throw OomphLibError(
442  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
443  break;
444  }
445 
446  // Dummy return (never get here)
447  return 0;
448  }
449 
450 
451  //==========================================================================
452  /// Return a combined error estimate from all compound flux errors
453  /// The default is to return the maximum of the compound flux errors
454  /// which will always force refinment if any field is above the single
455  /// mesh error threshold and unrefinement if both are below the lower limit.
456  /// Any other fancy combinations can be selected by
457  /// specifying a user-defined combined estimate by setting a function
458  /// pointer.
459  //==========================================================================
461  const Vector<double>& compound_error)
462  {
463  // If the function pointer has been set, call that function
464  if (Combined_error_fct_pt != 0)
465  {
466  return (*Combined_error_fct_pt)(compound_error);
467  }
468 
469  // Otherwise simply return the maximum of the compound errors
470  const unsigned n_compound_error = compound_error.size();
471 // If there are no errors then we have a problem
472 #ifdef PARANOID
473  if (n_compound_error == 0)
474  {
475  throw OomphLibError(
476  "No compound errors have been passed, so maximum cannot be found.",
477  OOMPH_CURRENT_FUNCTION,
478  OOMPH_EXCEPTION_LOCATION);
479  }
480 #endif
481 
482 
483  // Initialise the maxmimum to the first compound error
484  double max_error = compound_error[0];
485  // Loop over the other errors to find the maximum
486  //(we have already taken absolute values, so we don't need to do so here)
487  for (unsigned i = 1; i < n_compound_error; i++)
488  {
489  if (compound_error[i] > max_error)
490  {
491  max_error = compound_error[i];
492  }
493  }
494 
495  // Return the maximum
496  return max_error;
497  }
498 
499  //======================================================================
500  /// Setup patches: For each vertex node pointed to by nod_pt,
501  /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
502  /// contains the pointers to the elements that the node is part of.
503  /// Also returns a Vector of vertex nodes for use in get_element_errors.
504  //======================================================================
506  Mesh*& mesh_pt,
508  adjacent_elements_pt,
509  Vector<Node*>& vertex_node_pt)
510  {
511  // (see also note at the end of get_element_errors below)
512  // NOTE FOR FUTURE REFERENCE - revisit in case of adaptivity problems in
513  // parallel jobs at the boundaries between processes
514  //
515  // The current method for distributed problems neglects recovered flux
516  // contributions from patches that cannot be assembled from (vertex) nodes
517  // on the current process, but would be assembled if the problem was not
518  // distributed. The only nodes for which this is the case are nodes which
519  // lie on the exact boundary between processes.
520  //
521  // The suggested method for fixing this requires the current process (A) to
522  // receive information from the process (B) on which the patch is
523  // assembled. These patches are precisely the patches on process B which
524  // contain no halo elements. The contribution from these patches needs to
525  // be sent to the nodes (on process A) that are on the boundary between
526  // A and B. Therefore a map is required to denote such nodes; a node is on
527  // the boundary of a process if it is a member of at least one halo element
528  // and one non-halo element for that process. Create a vector of bools
529  // which is the size of the number of processes and make the entry true if
530  // the node (through a map(nod_pt,vector<bools> node_bnd)) is on the
531  // boundary for a process and false otherwise. This should be done here in
532  // setup_patches.
533  //
534  // When it comes to the error calculation in get_element_errors (see later)
535  // the communication needs to take place after rec_flux_map(nod_pt,i) has
536  // been assembled for all other patches. A separate vector for the patches
537  // to be sent needs to be assembled: the recovered flux contribution at
538  // a node on process B is sent to the equivalent node on process A if the
539  // patch contains ONLY halo elements AND the node is on the boundary between
540  // A and B (i.e. both the entries in the mapped vector are set to true).
541  //
542  // If anyone else read this and has any questions, I have a fuller
543  // explanation written out (with some relevant diagrams in the 2D case).
544  //
545  // Andrew.Gait@manchester.ac.uk
546 
547  // Auxiliary map that contains element-adjacency for ALL nodes
548  std::map<Node*, Vector<ElementWithZ2ErrorEstimator*>*>
549  aux_adjacent_elements_pt;
550 
551 #ifdef PARANOID
552  // Check if all elements request the same recovery order
553  unsigned ndisagree = 0;
554 #endif
555 
556  // Loop over all elements to setup adjacency for all nodes.
557  // Need to do this because midside nodes can be corner nodes for
558  // adjacent smaller elements! Admittedly, the inclusion of interior
559  // nodes is wasteful...
560  unsigned nelem = mesh_pt->nelement();
561  for (unsigned e = 0; e < nelem; e++)
562  {
564  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(e));
565 
566 #ifdef PARANOID
567  // Check if all elements request the same recovery order
568  if (el_pt->nrecovery_order() != Recovery_order)
569  {
570  ndisagree++;
571  }
572 #endif
573 
574  // Loop all nodes in element
575  unsigned nnod = el_pt->nnode();
576  for (unsigned n = 0; n < nnod; n++)
577  {
578  Node* nod_pt = el_pt->node_pt(n);
579 
580  // Has this node been considered before?
581  if (aux_adjacent_elements_pt[nod_pt] == 0)
582  {
583  // Create Vector of pointers to its adjacent elements
584  aux_adjacent_elements_pt[nod_pt] =
586  }
587 
588  // Add pointer to adjacent element
589  (*aux_adjacent_elements_pt[nod_pt]).push_back(el_pt);
590  // }
591  }
592  } // end element loop
593 
594 #ifdef PARANOID
595  // Check if all elements request the same recovery order
596  if (ndisagree != 0)
597  {
598  oomph_info
599  << "\n\n========================================================\n";
600  oomph_info << "WARNING: " << std::endl;
601  oomph_info << ndisagree << " out of " << mesh_pt->nelement()
602  << " elements\n";
603  oomph_info
604  << "have different preferences for the order of the recovery\n";
605  oomph_info << "shape functions. We are using: Recovery_order="
606  << Recovery_order << std::endl;
607  oomph_info
608  << "========================================================\n\n";
609  }
610 #endif
611 
612  // Loop over all elements, extract adjacency for corner nodes only
613  nelem = mesh_pt->nelement();
614  for (unsigned e = 0; e < nelem; e++)
615  {
617  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(e));
618 
619  // Loop over corner nodes
620  unsigned n_node = el_pt->nvertex_node();
621  for (unsigned n = 0; n < n_node; n++)
622  {
623  Node* nod_pt = el_pt->vertex_node_pt(n);
624 
625  // Has this node been considered before?
626  if (adjacent_elements_pt[nod_pt] == 0)
627  {
628  // Add the node pointer to the vertex node container
629  vertex_node_pt.push_back(nod_pt);
630 
631  // Create Vector of pointers to its adjacent elements
632  adjacent_elements_pt[nod_pt] =
634 
635  // Copy across:
636  unsigned nel = (*aux_adjacent_elements_pt[nod_pt]).size();
637  for (unsigned e = 0; e < nel; e++)
638  {
639  (*adjacent_elements_pt[nod_pt])
640  .push_back((*aux_adjacent_elements_pt[nod_pt])[e]);
641  }
642  }
643  }
644 
645  } // end of loop over elements
646 
647  // Cleanup
648  typedef std::map<Node*, Vector<ElementWithZ2ErrorEstimator*>*>::iterator
649  ITT;
650  for (ITT it = aux_adjacent_elements_pt.begin();
651  it != aux_adjacent_elements_pt.end();
652  it++)
653  {
654  delete it->second;
655  }
656  }
657 
658 
659  //======================================================================
660  /// Given the vector of elements that make up a patch,
661  /// the number of recovery and flux terms, and the
662  /// spatial dimension of the problem, compute
663  /// the matrix of recovered flux coefficients and return
664  /// a pointer to it.
665  //======================================================================
667  const Vector<ElementWithZ2ErrorEstimator*>& patch_el_pt,
668  const unsigned& num_recovery_terms,
669  const unsigned& num_flux_terms,
670  const unsigned& dim,
671  DenseMatrix<double>*& recovered_flux_coefficient_pt)
672  {
673  // Create/initialise matrix for linear system
674  DenseDoubleMatrix recovery_mat(num_recovery_terms, num_recovery_terms, 0.0);
675 
676  // Ceate/initialise vector of RHSs
677  Vector<Vector<double>> rhs(num_flux_terms);
678  for (unsigned irhs = 0; irhs < num_flux_terms; irhs++)
679  {
680  rhs[irhs].resize(num_recovery_terms);
681  for (unsigned j = 0; j < num_recovery_terms; j++)
682  {
683  rhs[irhs][j] = 0.0;
684  }
685  }
686 
687 
688  // Create a new integration scheme based on the recovery order
689  // in the elements
690  // Need to find the type of the element, default is to assume a quad
691  bool is_q_mesh = true;
692  // If we can dynamic cast to the TElementBase, then it's a triangle/tet
693  // Note that I'm assuming that all elements are of the same geometry, but
694  // if they weren't we could adapt...
695  if (dynamic_cast<TElementBase*>(patch_el_pt[0]))
696  {
697  is_q_mesh = false;
698  }
699 
700  Integral* const integ_pt = this->integral_rec(dim, is_q_mesh);
701 
702  // Loop over all elements in patch to assemble linear system
703  unsigned nelem = patch_el_pt.size();
704  for (unsigned e = 0; e < nelem; e++)
705  {
706  // Get pointer to element
707  ElementWithZ2ErrorEstimator* const el_pt = patch_el_pt[e];
708 
709  // Create storage for the recovery shape function values
710  Vector<double> psi_r(num_recovery_terms);
711 
712  // Create vector to hold local coordinates
713  Vector<double> s(dim);
714 
715  // Loop over the integration points
716  unsigned Nintpt = integ_pt->nweight();
717 
718  for (unsigned ipt = 0; ipt < Nintpt; ipt++)
719  {
720  // Assign values of s, the local coordinate
721  for (unsigned i = 0; i < dim; i++)
722  {
723  s[i] = integ_pt->knot(ipt, i);
724  }
725 
726  // Get the integral weight
727  double w = integ_pt->weight(ipt);
728 
729  // Jaocbian of mapping
730  double J = el_pt->J_eulerian(s);
731 
732  // Interpolate the global (Eulerian) coordinate
733  Vector<double> x(dim);
734  el_pt->interpolated_x(s, x);
735 
736 
737  // Premultiply the weights and the Jacobian
738  // and the geometric jacobian weight (used in axisymmetric
739  // and spherical coordinate systems)
740  double W = w * J * (el_pt->geometric_jacobian(x));
741 
742  // Recovery shape functions at global (Eulerian) coordinate
743  shape_rec(x, dim, psi_r);
744 
745  // Get FE estimates for Z2 flux:
746  Vector<double> fe_flux(num_flux_terms);
747  el_pt->get_Z2_flux(s, fe_flux);
748 
749  // Add elemental RHSs and recovery matrix to global versions
750  //----------------------------------------------------------
751 
752  // RHS for different flux components
753  for (unsigned i = 0; i < num_flux_terms; i++)
754  {
755  // Loop over the nodes for the test functions
756  for (unsigned l = 0; l < num_recovery_terms; l++)
757  {
758  rhs[i][l] += fe_flux[i] * psi_r[l] * W;
759  }
760  }
761 
762 
763  // Loop over the nodes for the test functions
764  for (unsigned l = 0; l < num_recovery_terms; l++)
765  {
766  // Loop over the nodes for the variables
767  for (unsigned l2 = 0; l2 < num_recovery_terms; l2++)
768  {
769  // Add contribution to recovery matrix
770  recovery_mat(l, l2) += psi_r[l] * psi_r[l2] * W;
771  }
772  }
773  }
774 
775  } // End of loop over elements that make up patch.
776 
777  // Delete the integration scheme
778  delete integ_pt;
779 
780  // Linear system is now assembled: Solve recovery system
781 
782  // LU decompose the recovery matrix
783  recovery_mat.ludecompose();
784 
785  // Back-substitute (and overwrite for all rhs
786  for (unsigned irhs = 0; irhs < num_flux_terms; irhs++)
787  {
788  recovery_mat.lubksub(rhs[irhs]);
789  }
790 
791  // Now create a matrix to store the flux recovery coefficients.
792  // Pointer to this matrix will be returned.
793  recovered_flux_coefficient_pt =
794  new DenseMatrix<double>(num_recovery_terms, num_flux_terms);
795 
796  // Copy coefficients
797  for (unsigned icoeff = 0; icoeff < num_recovery_terms; icoeff++)
798  {
799  for (unsigned irhs = 0; irhs < num_flux_terms; irhs++)
800  {
801  (*recovered_flux_coefficient_pt)(icoeff, irhs) = rhs[irhs][icoeff];
802  }
803  }
804  }
805 
806 
807  //==================================================================
808  /// Number of coefficients for expansion of recovered fluxes
809  /// for given spatial dimension of elements.
810  /// Use complete polynomial of given order for recovery
811  //==================================================================
812  unsigned Z2ErrorEstimator::nrecovery_terms(const unsigned& dim)
813  {
814  unsigned num_recovery_terms;
815 
816 #ifdef PARANOID
817  if ((dim != 1) && (dim != 2) && (dim != 3))
818  {
819  std::string error_message = "THIS HASN'T BEEN USED/VALIDATED YET -- "
820  "CHECK NUMBER OF RECOVERY TERMS!\n";
821  error_message += "Then remove this break and continue\n";
822 
823  throw OomphLibError(
824  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
825  }
826 #endif
827 
828  switch (Recovery_order)
829  {
830  case 1:
831 
832  // Linear recovery shape functions
833  //--------------------------------
834 
835  switch (dim)
836  {
837  case 1:
838  // 1D:
839  num_recovery_terms = 2; // 1, x
840  break;
841 
842  case 2:
843  // 2D:
844  num_recovery_terms = 3; // 1, x, y
845  break;
846 
847  case 3:
848  // 3D:
849  num_recovery_terms = 4; // 1, x, y, z
850  break;
851 
852  default:
853  throw OomphLibError("Dim must be 1, 2 or 3",
854  OOMPH_CURRENT_FUNCTION,
855  OOMPH_EXCEPTION_LOCATION);
856  }
857  break;
858 
859  case 2:
860 
861  // Quadratic recovery shape functions
862  //-----------------------------------
863 
864  switch (dim)
865  {
866  case 1:
867  // 1D:
868  num_recovery_terms = 3; // 1, x, x^2
869  break;
870 
871  case 2:
872  // 2D:
873  num_recovery_terms = 6; // 1, x, y, x^2, xy, y^2
874  break;
875 
876  case 3:
877  // 3D:
878  num_recovery_terms = 10; // 1, x, y, z, x^2, y^2, z^2, xy, xz, yz
879  break;
880 
881  default:
882  throw OomphLibError("Dim must be 1, 2 or 3",
883  OOMPH_CURRENT_FUNCTION,
884  OOMPH_EXCEPTION_LOCATION);
885  }
886  break;
887 
888 
889  case 3:
890 
891  // Cubic recovery shape functions
892  //--------------------------------
893 
894  switch (dim)
895  {
896  case 1:
897  // 1D:
898  num_recovery_terms = 4; // 1, x, x^2, x^3
899  break;
900 
901  case 2:
902  // 2D:
903  num_recovery_terms =
904  10; // 1, x, y, x^2, xy, y^2, x^3, y^3, x^2 y, x y^2
905  break;
906 
907  case 3:
908  // 3D:
909  num_recovery_terms = 20; // 1, x, y, z, x^2, y^2, z^2, xy, xz, yz,
910  // x^3, y^3, z^3,
911  // x^2 y, x^2 z,
912  // y^2 x, y^2 z,
913  // z^2 x, z^2 y
914  // xyz
915  break;
916 
917  default:
918  throw OomphLibError("Dim must be 1, 2 or 3",
919  OOMPH_CURRENT_FUNCTION,
920  OOMPH_EXCEPTION_LOCATION);
921  }
922  break;
923 
924 
925  default:
926 
927  // Any other recovery order?
928  //--------------------------
929  std::ostringstream error_stream;
930  error_stream << "Wrong Recovery_order " << Recovery_order << std::endl;
931 
932  throw OomphLibError(
933  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
934  }
935 
936  return num_recovery_terms;
937  }
938 
939 
940  //======================================================================
941  /// Get Vector of Z2-based error estimates for all elements in mesh.
942  /// If doc_info.is_doc_enabled()=true, doc FE and recovered fluxes in
943  /// - flux_fe*.dat
944  /// - flux_rec*.dat
945  //======================================================================
947  Vector<double>& elemental_error,
948  DocInfo& doc_info)
949  {
950 #ifdef OOMPH_HAS_MPI
951  // Storage for number of processors and current processor
952  int n_proc = 1;
953  int my_rank = 0;
954  if (mesh_pt->communicator_pt() != 0)
955  {
956  my_rank = mesh_pt->communicator_pt()->my_rank();
957  n_proc = mesh_pt->communicator_pt()->nproc();
958  }
960  {
961  my_rank = MPI_Helpers::communicator_pt()->my_rank();
962  n_proc = MPI_Helpers::communicator_pt()->nproc();
963  }
964 
965  MPI_Status status;
966 
967  // Initialise local values for all processes on mesh
968  unsigned num_flux_terms_local = 0;
969  unsigned dim_local = 0;
970  unsigned recovery_order_local = 0;
971 #endif
972 
973  // Global variables
974  unsigned num_flux_terms = 0;
975  unsigned dim = 0;
976 
977 #ifdef OOMPH_HAS_MPI
978  // It may be possible that a submesh contains no elements on a
979  // particular process after distribution. In order to instigate the
980  // error estimator calculations we need some information from the
981  // "first" element in a mesh; the following uses an MPI_Allreduce
982  // to figure out this information and communicate it to all processors
983  if (mesh_pt->nelement() > 0)
984  {
985  // Extract a few vital parameters from first element in mesh:
987  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(0));
988 #ifdef PARANOID
989  if (el_pt == 0)
990  {
991  throw OomphLibError(
992  "Element needs to inherit from ElementWithZ2ErrorEstimator!",
993  OOMPH_CURRENT_FUNCTION,
994  OOMPH_EXCEPTION_LOCATION);
995  }
996 #endif
997 
998  // Number of 'flux'-like terms to be recovered
999  num_flux_terms_local = el_pt->num_Z2_flux_terms();
1000 
1001  // Determine spatial dimension of all elements from first element
1002  dim_local = el_pt->dim();
1003 
1004  // Do we need to determine the recovery order from first element?
1006  {
1007  recovery_order_local = el_pt->nrecovery_order();
1008  }
1009 
1010  } // end if (mesh_pt->nelement()>0)
1011 
1012  // Storage for the recovery order
1013  unsigned recovery_order = 0;
1014 
1015  // Now communicate these via an MPI_Allreduce to every process
1016  // if the mesh has been distributed
1017  if (mesh_pt->is_mesh_distributed())
1018  {
1019  // Get communicator from mesh
1020  OomphCommunicator* comm_pt = mesh_pt->communicator_pt();
1021 
1022  MPI_Allreduce(&num_flux_terms_local,
1023  &num_flux_terms,
1024  1,
1025  MPI_UNSIGNED,
1026  MPI_MAX,
1027  comm_pt->mpi_comm());
1028  MPI_Allreduce(&dim_local, &dim, 1, MPI_INT, MPI_MAX, comm_pt->mpi_comm());
1029  MPI_Allreduce(&recovery_order_local,
1030  &recovery_order,
1031  1,
1032  MPI_UNSIGNED,
1033  MPI_MAX,
1034  comm_pt->mpi_comm());
1035  }
1036  else
1037  {
1038  num_flux_terms = num_flux_terms_local;
1039  dim = dim_local;
1040  recovery_order = recovery_order_local;
1041  }
1042 
1043  // Do we need to determine the recovery order from first element?
1045  {
1047  }
1048 
1049 #else // !OOMPH_HAS_MPI
1050 
1051  // Extract a few vital parameters from first element in mesh:
1053  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(0));
1054 #ifdef PARANOID
1055  if (el_pt == 0)
1056  {
1057  throw OomphLibError(
1058  "Element needs to inherit from ElementWithZ2ErrorEstimator!",
1059  OOMPH_CURRENT_FUNCTION,
1060  OOMPH_EXCEPTION_LOCATION);
1061  }
1062 #endif
1063 
1064  // Number of 'flux'-like terms to be recovered
1065  num_flux_terms = el_pt->num_Z2_flux_terms();
1066 
1067  // Determine spatial dimension of all elements from first element
1068  dim = el_pt->dim();
1069 
1070  // Do we need to determine the recovery order from first element?
1072  {
1073  Recovery_order = el_pt->nrecovery_order();
1074  }
1075 
1076 #endif
1077 
1078  // Determine number of coefficients for expansion of recovered fluxes
1079  // Use complete polynomial of given order for recovery
1080  unsigned num_recovery_terms = nrecovery_terms(dim);
1081 
1082 
1083  // Setup patches (also returns Vector of vertex nodes)
1084  //====================================================
1085  std::map<Node*, Vector<ElementWithZ2ErrorEstimator*>*> adjacent_elements_pt;
1086  Vector<Node*> vertex_node_pt;
1087  setup_patches(mesh_pt, adjacent_elements_pt, vertex_node_pt);
1088 
1089  // Loop over all patches to get recovered flux value coefficients
1090  //===============================================================
1091 
1092  // Map to store sets of pointers to the recovered flux coefficient matrices
1093  // for each node.
1094  std::map<Node*, std::set<DenseMatrix<double>*>> flux_coeff_pt;
1095 
1096  // We store the pointers to the recovered flux coefficient matrices for
1097  // various patches in a vector so we can delete them later
1098  Vector<DenseMatrix<double>*> vector_of_recovered_flux_coefficient_pt;
1099 
1100  typedef std::map<Node*, Vector<ElementWithZ2ErrorEstimator*>*>::iterator IT;
1101 
1102  // Need to translate ElementWithZ2ErrorEstimator pointer to element number
1103  // in order to give each processor elements to work on if the problem
1104  // has not yet been distributed. In order to reduce the use of #ifdef this
1105  // is also done for the serial problem and the code is amended accordingly.
1106 
1107  std::map<ElementWithZ2ErrorEstimator*, int> elem_num;
1108  unsigned nelem = mesh_pt->nelement();
1109  for (unsigned e = 0; e < nelem; e++)
1110  {
1111  elem_num[dynamic_cast<ElementWithZ2ErrorEstimator*>(
1112  mesh_pt->element_pt(e))] = e;
1113  }
1114 
1115  // This isn't a global variable
1116  int n_patch = adjacent_elements_pt.size(); // also needed by serial version
1117 
1118  // Default values for serial AND parallel distributed problem
1119  int itbegin = 0;
1120 
1121  int itend = n_patch;
1122 
1123 #ifdef OOMPH_HAS_MPI
1124  int range = n_patch;
1125 
1126  // Work out values for parallel non-distributed problem
1127  if (!(mesh_pt->is_mesh_distributed()))
1128  {
1129  // setup the loop variables
1130  range = n_patch / n_proc; // number of patches on each proc
1131 
1132  itbegin = my_rank * range;
1133  itend = (my_rank + 1) * range;
1134 
1135  // if on the last processor, ensure the end matches
1136  if (my_rank == (n_proc - 1))
1137  {
1138  itend = n_patch;
1139  }
1140  }
1141 #endif
1142 
1143  // Set up matrices and vectors which will be sent later
1144  // - full matrix of all recovered coefficients
1146  vector_of_recovered_flux_coefficient_pt_to_send;
1147  // - vectors containing element numbers in each patch
1148  Vector<Vector<int>> vector_of_elements_in_patch_to_send;
1149 
1150  // Now we can loop over the patches on the current process
1151  for (int i = itbegin; i < itend; i++)
1152  {
1153  // Which vertex node are we at?
1154  Node* nod_pt = vertex_node_pt[i];
1155 
1156  // Pointer to vector of pointers to elements that make up
1157  // the patch.
1159  adjacent_elements_pt[nod_pt];
1160 
1161  // Is the corner node that is central to the patch surrounded by
1162  // at least two elements?
1163  unsigned nelem = (*el_vec_pt).size();
1164 
1165  if (nelem >= 2)
1166  {
1167  // store the number of elements in the patch
1168  Vector<int> elements_in_this_patch;
1169  for (unsigned e = 0; e < nelem; e++)
1170  {
1171  elements_in_this_patch.push_back(elem_num[(*el_vec_pt)[e]]);
1172  }
1173 
1174  // put them into storage vector ready to send
1175  vector_of_elements_in_patch_to_send.push_back(elements_in_this_patch);
1176 
1177  // Given the vector of elements that make up the patch,
1178  // the number of recovery and flux terms, and the spatial
1179  // dimension of the problem, compute
1180  // the matrix of recovered flux coefficients and return
1181  // a pointer to it.
1182  DenseMatrix<double>* recovered_flux_coefficient_pt = 0;
1183  get_recovered_flux_in_patch(*el_vec_pt,
1184  num_recovery_terms,
1185  num_flux_terms,
1186  dim,
1187  recovered_flux_coefficient_pt);
1188 
1189  // Store pointer to recovered flux coefficients for
1190  // current patch in vector so we can send and then delete it later
1191  vector_of_recovered_flux_coefficient_pt_to_send.push_back(
1192  recovered_flux_coefficient_pt);
1193 
1194  } // end of if(nelem>=2)
1195  }
1196 
1197  // Now broadcast the result from each process to every other process
1198  // if the mesh has not yet been distributed and MPI is initialised
1199 
1200 #ifdef OOMPH_HAS_MPI
1201 
1202  if (!mesh_pt->is_mesh_distributed() &&
1204  {
1205  // Get communicator from namespace
1207 
1208  // All local recovered fluxes have been calculated, so now share result
1209  for (int iproc = 0; iproc < n_proc; iproc++)
1210  {
1211  // Broadcast number of patches processed
1212  int n_patches = vector_of_recovered_flux_coefficient_pt_to_send.size();
1213  MPI_Bcast(&n_patches, 1, MPI_INT, iproc, comm_pt->mpi_comm());
1214 
1215  // Loop over these patches, broadcast recovered flux coefficients
1216  for (int ipatch = 0; ipatch < n_patches; ipatch++)
1217  {
1218  // Number of elements in this patch
1219  Vector<int> elements(0);
1220  unsigned nelements = 0;
1221 
1222  // Which processor are we on?
1223  if (my_rank == iproc)
1224  {
1225  elements = vector_of_elements_in_patch_to_send[ipatch];
1226  nelements = elements.size();
1227  }
1228 
1229  // Broadcast elements
1230  comm_pt->broadcast(iproc, elements);
1231 
1232  // Now get recovered flux coefficients
1233  DenseMatrix<double>* recovered_flux_coefficient_pt;
1234 
1235  // Which processor are we on?
1236  if (my_rank == iproc)
1237  {
1238  recovered_flux_coefficient_pt =
1239  vector_of_recovered_flux_coefficient_pt_to_send[ipatch];
1240  }
1241  else
1242  {
1243  recovered_flux_coefficient_pt = new DenseMatrix<double>;
1244  }
1245 
1246  // broadcast this matrix from the loop processor
1247  DenseMatrix<double> mattosend = *recovered_flux_coefficient_pt;
1248  comm_pt->broadcast(iproc, mattosend);
1249 
1250  // Set pointer on all processors
1251  *recovered_flux_coefficient_pt = mattosend;
1252 
1253  // End of parallel broadcasting
1254  vector_of_recovered_flux_coefficient_pt.push_back(
1255  recovered_flux_coefficient_pt);
1256 
1257  // Loop over elements in patch (work out nelements again after bcast)
1258  nelements = elements.size();
1259 
1260  for (unsigned e = 0; e < nelements; e++)
1261  {
1262  // Get pointer to element
1264  dynamic_cast<ElementWithZ2ErrorEstimator*>(
1265  mesh_pt->element_pt(elements[e]));
1266 
1267  // Loop over all nodes in element
1268  unsigned num_nod = el_pt->nnode();
1269  for (unsigned n = 0; n < num_nod; n++)
1270  {
1271  // Get the node
1272  Node* nod_pt = el_pt->node_pt(n);
1273  // Add the pointer to the current flux coefficient matrix
1274  // to the set for the node
1275  // Mesh not distributed here so nod_pt cannot be halo
1276  flux_coeff_pt[nod_pt].insert(recovered_flux_coefficient_pt);
1277  }
1278  }
1279 
1280  } // end loop over patches on current processor
1281 
1282  } // end loop over processors
1283  }
1284  else // is_mesh_distributed=true
1285  {
1286 #endif // end ifdef OOMPH_HAS_MPI for parallel job without mesh distribution
1287 
1288  // Do the same for a distributed mesh as for a serial job
1289  // up to the point where the elemental error is calculated
1290  // and then communicate that (see below)
1291  int n_patches = vector_of_recovered_flux_coefficient_pt_to_send.size();
1292 
1293  // Loop over these patches
1294  for (int ipatch = 0; ipatch < n_patches; ipatch++)
1295  {
1296  // Number of elements in this patch
1297  Vector<int> elements;
1298  int nelements; // Needs to be int for elem_num call later
1299  elements = vector_of_elements_in_patch_to_send[ipatch];
1300  nelements = elements.size();
1301 
1302  // Now get recovered flux coefficients
1303  DenseMatrix<double>* recovered_flux_coefficient_pt;
1304  recovered_flux_coefficient_pt =
1305  vector_of_recovered_flux_coefficient_pt_to_send[ipatch];
1306 
1307  vector_of_recovered_flux_coefficient_pt.push_back(
1308  recovered_flux_coefficient_pt);
1309 
1310  for (int e = 0; e < nelements; e++)
1311  {
1312  // Get pointer to element
1314  dynamic_cast<ElementWithZ2ErrorEstimator*>(
1315  mesh_pt->element_pt(elements[e]));
1316 
1317  // Loop over all nodes in element
1318  unsigned num_nod = el_pt->nnode();
1319  for (unsigned n = 0; n < num_nod; n++)
1320  {
1321  // Get the node
1322  Node* nod_pt = el_pt->node_pt(n);
1323  // Add the pointer to the current flux coefficient matrix
1324  // to the set for this node
1325  flux_coeff_pt[nod_pt].insert(recovered_flux_coefficient_pt);
1326  }
1327  }
1328 
1329  } // End loop over patches on current processor
1330 
1331 
1332 #ifdef OOMPH_HAS_MPI
1333  } // End if(is_mesh_distributed)
1334 #endif
1335 
1336  // Cleanup patch storage scheme
1337  for (IT it = adjacent_elements_pt.begin(); it != adjacent_elements_pt.end();
1338  it++)
1339  {
1340  delete it->second;
1341  }
1342  adjacent_elements_pt.clear();
1343 
1344  // Loop over all nodes, take average of recovered flux values
1345  //-----------------------------------------------------------
1346  // and evaluate recovered flux at nodes
1347  //-------------------------------------
1348 
1349  // Map of (averaged) recoverd flux values at nodes
1351 
1352  // Loop over all nodes
1353  unsigned n_node = mesh_pt->nnode();
1354  for (unsigned n = 0; n < n_node; n++)
1355  {
1356  Node* nod_pt = mesh_pt->node_pt(n);
1357 
1358  // How many patches is this node a member of?
1359  unsigned npatches = flux_coeff_pt[nod_pt].size();
1360 
1361  // Matrix of averaged coefficients for this node
1362  DenseMatrix<double> averaged_flux_coeff(
1363  num_recovery_terms, num_flux_terms, 0.0);
1364 
1365  // Loop over matrices for different patches and add contributions
1366  typedef std::set<DenseMatrix<double>*>::iterator IT;
1367  for (IT it = flux_coeff_pt[nod_pt].begin();
1368  it != flux_coeff_pt[nod_pt].end();
1369  it++)
1370  {
1371  for (unsigned i = 0; i < num_recovery_terms; i++)
1372  {
1373  for (unsigned j = 0; j < num_flux_terms; j++)
1374 
1375  {
1376  // ...just add it -- we'll divide by the number of patches later
1377  averaged_flux_coeff(i, j) += (*(*it))(i, j);
1378  }
1379  }
1380  }
1381 
1382  // Now evaluate the recovered flux (based on the averaged coefficients)
1383  //---------------------------------------------------------------------
1384  // at the nodal position itself.
1385  //------------------------------
1386 
1387  // Get global (Eulerian) nodal position
1388  Vector<double> x(dim);
1389  for (unsigned i = 0; i < dim; i++)
1390  {
1391  x[i] = nod_pt->x(i);
1392  }
1393 
1394  // Evaluate global recovery functions at node
1395  Vector<double> psi_r(num_recovery_terms);
1396  shape_rec(x, dim, psi_r);
1397 
1398  // Initialise nodal fluxes
1399  for (unsigned i = 0; i < num_flux_terms; i++)
1400  {
1401  rec_flux_map(nod_pt, i) = 0.0;
1402  }
1403 
1404  // Loop over coefficients for flux recovery
1405  for (unsigned i = 0; i < num_flux_terms; i++)
1406  {
1407  for (unsigned icoeff = 0; icoeff < num_recovery_terms; icoeff++)
1408  {
1409  rec_flux_map(nod_pt, i) +=
1410  averaged_flux_coeff(icoeff, i) * psi_r[icoeff];
1411  }
1412  // Now take averaging into account
1413  rec_flux_map(nod_pt, i) /= double(npatches);
1414  }
1415 
1416  } // end loop over nodes
1417 
1418  // We're done with the recovered flux coefficient matrices for
1419  // the various patches and can delete them
1420  unsigned npatch = vector_of_recovered_flux_coefficient_pt.size();
1421  for (unsigned p = 0; p < npatch; p++)
1422  {
1423  delete vector_of_recovered_flux_coefficient_pt[p];
1424  }
1425 
1426  // NOTE FOR FUTURE REFERENCE - revisit in case of adaptivity problems in
1427  // parallel jobs
1428  //
1429  // The current method for distributed problems neglects recovered flux
1430  // contributions from patches that cannot be assembled on the current
1431  // processor, but would be assembled if the problem was not distributed.
1432  // The only nodes for which this is the case are nodes which lie on the
1433  // exact boundary between processors.
1434  //
1435  // See the note at the start of setup_patches (above) for more details.
1436 
1437  // Get error estimates for all elements
1438  //======================================
1439 
1440  // Find the number of compound fluxes
1441  // Loop over all (non-halo) elements
1442  nelem = mesh_pt->nelement();
1443  // Initialise the number of compound fluxes
1444  // Must be an integer for an MPI call later on
1445  int n_compound_flux = 1;
1446  for (unsigned e = 0; e < nelem; e++)
1447  {
1449  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(e));
1450 
1451 #ifdef OOMPH_HAS_MPI
1452  // Ignore halo elements
1453  if (!el_pt->is_halo())
1454  {
1455 #endif
1456  // Find the number of compound fluxes in the element
1457  const int n_compound_flux_el = el_pt->ncompound_fluxes();
1458  // If it's greater than the current (global) number of compound fluxes
1459  // bump up the global number
1460  if (n_compound_flux_el > n_compound_flux)
1461  {
1462  n_compound_flux = n_compound_flux_el;
1463  }
1464 #ifdef OOMPH_HAS_MPI
1465  } // end if (!el_pt->is_halo())
1466 #endif
1467  }
1468 
1469  // Initialise a vector of flux norms
1470  Vector<double> flux_norm(n_compound_flux, 0.0);
1471 
1472  unsigned test_count = 0;
1473 
1474  // Storage for the elemental compound flux error
1475  DenseMatrix<double> elemental_compound_flux_error(
1476  nelem, n_compound_flux, 0.0);
1477 
1478  // Loop over all (non-halo) elements again
1479  for (unsigned e = 0; e < nelem; e++)
1480  {
1482  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(e));
1483 
1484 #ifdef OOMPH_HAS_MPI
1485  // Ignore halo elements
1486  if (!el_pt->is_halo())
1487  {
1488 #endif
1489 
1490  Vector<double> s(dim);
1491 
1492  // Initialise elemental error one for each compound flux in the element
1493  const unsigned n_compound_flux_el = el_pt->ncompound_fluxes();
1494  Vector<double> error(n_compound_flux_el, 0.0);
1495 
1496  Integral* integ_pt = el_pt->integral_pt();
1497 
1498  // Set the value of Nintpt
1499  const unsigned n_intpt = integ_pt->nweight();
1500 
1501  // Loop over the integration points
1502  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1503  {
1504  // Assign values of s
1505  for (unsigned i = 0; i < dim; i++)
1506  {
1507  s[i] = integ_pt->knot(ipt, i);
1508  }
1509 
1510  // Get the integral weight
1511  double w = integ_pt->weight(ipt);
1512 
1513  // Jacobian of mapping
1514  double J = el_pt->J_eulerian(s);
1515 
1516  // Get the Eulerian position
1517  Vector<double> x(dim);
1518  el_pt->interpolated_x(s, x);
1519 
1520  // Premultiply the weights and the Jacobian
1521  // and the geometric jacobian weight (used in axisymmetric
1522  // and spherical coordinate systems)
1523  double W = w * J * (el_pt->geometric_jacobian(x));
1524 
1525  // Number of FE nodes
1526  unsigned n_node = el_pt->nnode();
1527 
1528  // FE shape function
1529  Shape psi(n_node);
1530 
1531  // Get values of FE shape function
1532  el_pt->shape(s, psi);
1533 
1534  // Initialise recovered flux Vector
1535  Vector<double> rec_flux(num_flux_terms, 0.0);
1536 
1537  // Loop over all nodes (incl. halo nodes) to assemble contribution
1538  for (unsigned n = 0; n < n_node; n++)
1539  {
1540  Node* nod_pt = el_pt->node_pt(n);
1541 
1542  // Loop over components
1543  for (unsigned i = 0; i < num_flux_terms; i++)
1544  {
1545  rec_flux[i] += rec_flux_map(nod_pt, i) * psi[n];
1546  }
1547  }
1548 
1549  // FE flux
1550  Vector<double> fe_flux(num_flux_terms);
1551  el_pt->get_Z2_flux(s, fe_flux);
1552  // Get compound flux indices. Initialised to zero
1553  Vector<unsigned> flux_index(num_flux_terms, 0);
1554  el_pt->get_Z2_compound_flux_indices(flux_index);
1555 
1556  // Add to RMS errors for each compound flux:
1557  Vector<double> sum(n_compound_flux_el, 0.0);
1558  Vector<double> sum2(n_compound_flux_el, 0.0);
1559  for (unsigned i = 0; i < num_flux_terms; i++)
1560  {
1561  sum[flux_index[i]] +=
1562  (rec_flux[i] - fe_flux[i]) * (rec_flux[i] - fe_flux[i]);
1563  sum2[flux_index[i]] += rec_flux[i] * rec_flux[i];
1564  }
1565 
1566  for (unsigned i = 0; i < n_compound_flux_el; i++)
1567  {
1568  // Add the errors to the appropriate compound flux error
1569  error[i] += sum[i] * W;
1570  // Add to flux norm
1571  flux_norm[i] += sum2[i] * W;
1572  }
1573  }
1574  // Unscaled elemental RMS error:
1575  test_count++; // counting elements visited
1576 
1577  // elemental_error[e]=sqrt(error);
1578  // Take the square-root of the appropriate flux error and
1579  // store the result
1580  for (unsigned i = 0; i < n_compound_flux_el; i++)
1581  {
1582  elemental_compound_flux_error(e, i) = sqrt(error[i]);
1583  }
1584 
1585 #ifdef OOMPH_HAS_MPI
1586  } // end if (!el_pt->is_halo())
1587 #endif
1588  } // end of loop over elements
1589 
1590  // Communicate the error for haloed elements to halo elements:
1591  // - loop over processors
1592  // - if current process, receive to halo element error
1593  // - if not current process, send haloed element error
1594  // How do we know which part of elemental_error to send?
1595  // Loop over haloed elements and find element number using elem_num map;
1596  // send this - order preservation of halo/haloed elements
1597  // guarantees that they get through in the correct order
1598 #ifdef OOMPH_HAS_MPI
1599 
1600  if (mesh_pt->is_mesh_distributed())
1601  {
1602  // Get communicator from mesh
1603  OomphCommunicator* comm_pt = mesh_pt->communicator_pt();
1604 
1605  for (int iproc = 0; iproc < n_proc; iproc++)
1606  {
1607  if (iproc != my_rank) // Not current process, so send
1608  {
1609  // Get the haloed elements
1610  Vector<GeneralisedElement*> haloed_elem_pt =
1611  mesh_pt->haloed_element_pt(iproc);
1612  // Find the number of haloed elements
1613  int nelem_haloed = haloed_elem_pt.size();
1614 
1615  // If there are some haloed elements, assemble and send the
1616  // errors
1617  if (nelem_haloed != 0)
1618  {
1619  // Find the number of error entires:
1620  // number of haloed elements x number of compound fluxes
1621  int n_elem_error_haloed = nelem_haloed * n_compound_flux;
1622  // Vector for elemental errors
1623  Vector<double> haloed_elem_error(n_elem_error_haloed);
1624  // Counter for the vector index
1625  unsigned count = 0;
1626  for (int e = 0; e < nelem_haloed; e++)
1627  {
1628  // Find element number
1629  int element_num =
1630  elem_num[dynamic_cast<ElementWithZ2ErrorEstimator*>(
1631  haloed_elem_pt[e])];
1632  // Put the error in a vector to send
1633  for (int i = 0; i < n_compound_flux; i++)
1634  {
1635  haloed_elem_error[count] =
1636  elemental_compound_flux_error(element_num, i);
1637  ++count;
1638  }
1639  }
1640  // Send the errors
1641  MPI_Send(&haloed_elem_error[0],
1642  n_elem_error_haloed,
1643  MPI_DOUBLE,
1644  iproc,
1645  0,
1646  comm_pt->mpi_comm());
1647  }
1648  }
1649  else // iproc=my_rank, so receive errors from others
1650  {
1651  for (int send_rank = 0; send_rank < n_proc; send_rank++)
1652  {
1653  if (iproc != send_rank) // iproc=my_rank already!
1654  {
1655  Vector<GeneralisedElement*> halo_elem_pt =
1656  mesh_pt->halo_element_pt(send_rank);
1657  // Find number of halo elements
1658  int nelem_halo = halo_elem_pt.size();
1659  // If there are some halo elements, receive errors and
1660  // put in the appropriate places
1661  if (nelem_halo != 0)
1662  {
1663  // Find the number of error entires:
1664  // number of haloed elements x number of compound fluxes
1665  int n_elem_error_halo = nelem_halo * n_compound_flux;
1666  // Vector for elemental errors
1667  Vector<double> halo_elem_error(n_elem_error_halo);
1668 
1669 
1670  // Receive the errors from processor send_rank
1671  MPI_Recv(&halo_elem_error[0],
1672  n_elem_error_halo,
1673  MPI_DOUBLE,
1674  send_rank,
1675  0,
1676  comm_pt->mpi_comm(),
1677  &status);
1678 
1679  // Counter for the vector index
1680  unsigned count = 0;
1681  for (int e = 0; e < nelem_halo; e++)
1682  {
1683  // Find element number
1684  int element_num =
1685  elem_num[dynamic_cast<ElementWithZ2ErrorEstimator*>(
1686  halo_elem_pt[e])];
1687  // Put the error in the correct location
1688  for (int i = 0; i < n_compound_flux; i++)
1689  {
1690  elemental_compound_flux_error(element_num, i) =
1691  halo_elem_error[count];
1692  ++count;
1693  }
1694  }
1695  }
1696  }
1697  } // End of interior loop over processors
1698  }
1699  } // End of exterior loop over processors
1700 
1701  } // End of if (mesh has been distributed)
1702 
1703 #endif
1704 
1705  // NOTE FOR FUTURE REFERENCE - revisit in case of adaptivity problems in
1706  // parallel jobs
1707  //
1708  // The current method for distributed problems neglects recovered flux
1709  // contributions from patches that cannot be assembled on the current
1710  // processor, but would be assembled if the problem was not distributed.
1711  // The only nodes for which this is the case are nodes which lie on the
1712  // exact boundary between processors.
1713  //
1714  // See the note at the start of setup_patches (above) for more details.
1715 
1716  // Use computed flux norm or externally imposed reference value?
1717  if (Reference_flux_norm != 0.0)
1718  {
1719  // At the moment assume that all fluxes have the same reference norm
1720  for (int i = 0; i < n_compound_flux; i++)
1721  {
1722  flux_norm[i] = Reference_flux_norm;
1723  }
1724  }
1725  else
1726  {
1727  // In parallel, perform reduction operation to get global value
1728 #ifdef OOMPH_HAS_MPI
1729  if (mesh_pt->is_mesh_distributed())
1730  {
1731  // Get communicator from mesh
1732  OomphCommunicator* comm_pt = mesh_pt->communicator_pt();
1733 
1734  Vector<double> total_flux_norm(n_compound_flux);
1735  // every process needs to know the sum
1736  MPI_Allreduce(&flux_norm[0],
1737  &total_flux_norm[0],
1738  n_compound_flux,
1739  MPI_DOUBLE,
1740  MPI_SUM,
1741  comm_pt->mpi_comm());
1742  // take sqrt
1743  for (int i = 0; i < n_compound_flux; i++)
1744  {
1745  flux_norm[i] = sqrt(total_flux_norm[i]);
1746  }
1747  }
1748  else // mesh has not been distributed, so flux_norm already global
1749  {
1750  for (int i = 0; i < n_compound_flux; i++)
1751  {
1752  flux_norm[i] = sqrt(flux_norm[i]);
1753  }
1754  }
1755 #else // serial problem, so flux_norm already global
1756  for (int i = 0; i < n_compound_flux; i++)
1757  {
1758  flux_norm[i] = sqrt(flux_norm[i]);
1759  }
1760 #endif
1761  }
1762 
1763  // Now loop over (all!) elements again and
1764  // normalise errors by global flux norm
1765 
1766  nelem = mesh_pt->nelement();
1767  for (unsigned e = 0; e < nelem; e++)
1768  {
1769  // Get the vector of normalised compound fluxes
1770  Vector<double> normalised_compound_flux_error(n_compound_flux);
1771  for (int i = 0; i < n_compound_flux; i++)
1772  {
1773  if (flux_norm[i] != 0.0)
1774  {
1775  normalised_compound_flux_error[i] =
1776  elemental_compound_flux_error(e, i) / flux_norm[i];
1777  }
1778  else
1779  {
1780  normalised_compound_flux_error[i] =
1781  elemental_compound_flux_error(e, i);
1782  }
1783  }
1784 
1785  // calculate the combined error estimate
1786  elemental_error[e] =
1787  get_combined_error_estimate(normalised_compound_flux_error);
1788  }
1789 
1790  // Doc global fluxes?
1791  if (doc_info.is_doc_enabled())
1792  {
1793  doc_flux(
1794  mesh_pt, num_flux_terms, rec_flux_map, elemental_error, doc_info);
1795  }
1796  }
1797 
1798 
1799  //==================================================================
1800  /// Doc FE and recovered flux
1801  //==================================================================
1803  Mesh* mesh_pt,
1804  const unsigned& num_flux_terms,
1805  MapMatrixMixed<Node*, int, double>& rec_flux_map,
1806  const Vector<double>& elemental_error,
1807  DocInfo& doc_info)
1808  {
1809 #ifdef OOMPH_HAS_MPI
1810 
1811  // Get communicator from mesh
1812  OomphCommunicator* comm_pt = mesh_pt->communicator_pt();
1813 
1814 #else
1815 
1816  // Dummy communicator
1818 
1819 #endif
1820 
1821  // File suffix identifying processor rank. If comm_pt is null (because
1822  // oomph-lib was built with MPI but this mesh is not distributed) the
1823  // string is empty.
1824  std::string rank_string = "";
1825  if (comm_pt != 0)
1826  {
1827  rank_string = "_on_proc_" + comm_pt->my_rank();
1828  }
1829 
1830  // Setup output files
1831  std::ofstream some_file, feflux_file;
1832  std::ostringstream filename;
1833  filename << doc_info.directory() << "/flux_rec" << doc_info.number()
1834  << rank_string << ".dat";
1835  some_file.open(filename.str().c_str());
1836  filename.str("");
1837  filename << doc_info.directory() << "/flux_fe" << doc_info.number()
1838  << rank_string << ".dat";
1839  feflux_file.open(filename.str().c_str());
1840 
1841  unsigned nel = mesh_pt->nelement();
1842  if (nel > 0)
1843  {
1844  // Extract first element to determine spatial dimension
1845  FiniteElement* el_pt = mesh_pt->finite_element_pt(0);
1846  unsigned dim = el_pt->dim();
1847  Vector<double> s(dim);
1848 
1849  // Decide on the number of plot points
1850  unsigned nplot = 5;
1851 
1852  // Loop over all elements
1853  for (unsigned e = 0; e < nel; e++)
1854  {
1856  dynamic_cast<ElementWithZ2ErrorEstimator*>(mesh_pt->element_pt(e));
1857 
1858  // Write tecplot header
1859  feflux_file << el_pt->tecplot_zone_string(nplot);
1860  some_file << el_pt->tecplot_zone_string(nplot);
1861 
1862  unsigned num_plot_points = el_pt->nplot_points(nplot);
1863  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1864  {
1865  // Get local coordinates of plot point
1866  el_pt->get_s_plot(iplot, nplot, s);
1867 
1868  // Coordinate
1869  Vector<double> x(dim);
1870  el_pt->interpolated_x(s, x);
1871 
1872  // Number of FE nodes
1873  unsigned n_node = el_pt->nnode();
1874 
1875  // FE shape function
1876  Shape psi(n_node);
1877 
1878  // Get values of FE shape function
1879  el_pt->shape(s, psi);
1880 
1881  // Initialise recovered flux Vector
1882  Vector<double> rec_flux(num_flux_terms, 0.0);
1883 
1884  // Loop over nodes to assemble contribution
1885  for (unsigned n = 0; n < n_node; n++)
1886  {
1887  Node* nod_pt = el_pt->node_pt(n);
1888 
1889  // Loop over components
1890  for (unsigned i = 0; i < num_flux_terms; i++)
1891  {
1892  rec_flux[i] += rec_flux_map(nod_pt, i) * psi[n];
1893  }
1894  }
1895 
1896  // FE flux
1897  Vector<double> fe_flux(num_flux_terms);
1898  el_pt->get_Z2_flux(s, fe_flux);
1899 
1900  for (unsigned i = 0; i < dim; i++)
1901  {
1902  some_file << x[i] << " ";
1903  }
1904  for (unsigned i = 0; i < num_flux_terms; i++)
1905  {
1906  some_file << rec_flux[i] << " ";
1907  }
1908  some_file << elemental_error[e] << " " << std::endl;
1909 
1910 
1911  for (unsigned i = 0; i < dim; i++)
1912  {
1913  feflux_file << x[i] << " ";
1914  }
1915  for (unsigned i = 0; i < num_flux_terms; i++)
1916  {
1917  feflux_file << fe_flux[i] << " ";
1918  }
1919  feflux_file << elemental_error[e] << " " << std::endl;
1920  }
1921 
1922  // Write tecplot footer (e.g. FE connectivity)
1923  // Do this for each element so the output is compatible with
1924  // oomph-convert
1925  el_pt->write_tecplot_zone_footer(some_file, nplot);
1926  el_pt->write_tecplot_zone_footer(feflux_file, nplot);
1927  }
1928  }
1929 
1930  some_file.close();
1931  feflux_file.close();
1932  }
1933 
1934 
1935 } // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
virtual void lubksub(DoubleVector &rhs)
LU backsubstitution.
Definition: matrices.cc:202
virtual void ludecompose()
LU decomposition using DenseLU (default linea solver)
Definition: matrices.cc:192
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
virtual unsigned ncompound_fluxes()
A stuitable error estimator for a multi-physics elements may require one Z2 error estimate for each f...
virtual void get_Z2_compound_flux_indices(Vector< unsigned > &flux_index)
Return the compound flux index of each flux component The default (do nothing behaviour) will mean th...
virtual void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)=0
Z2 'flux' terms for Z2 error estimation.
virtual double geometric_jacobian(const Vector< double > &x)
Return the geometric jacobian (should be overloaded in cylindrical and spherical geometries)....
virtual unsigned num_Z2_flux_terms()=0
Number of 'flux' terms for Z2 error estimation.
virtual unsigned nrecovery_order()=0
Order of recovery shape functions.
A general Finite Element class.
Definition: elements.h:1317
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4133
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3165
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2495
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3992
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2615
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1967
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3152
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2504
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3190
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3178
1D Gaussian integration class. Two integration points. This integration scheme can integrate up to th...
Definition: integral.h:159
1D Gaussian integration class. Three integration points. This integration scheme can integrate up to ...
Definition: integral.h:205
1D Gaussian integration class Four integration points. This integration scheme can integrate up to se...
Definition: integral.h:251
2D Gaussian integration class. 2x2 integration points. This integration scheme can integrate up to th...
Definition: integral.h:297
2D Gaussian integration class. 3x3 integration points. This integration scheme can integrate up to fi...
Definition: integral.h:343
2D Gaussian integration class. 4x4 integration points. This integration scheme can integrate up to se...
Definition: integral.h:389
3D Gaussian integration class 2x2x2 integration points. This integration scheme can integrate up to t...
Definition: integral.h:436
3D Gaussian integration class 3x3x3 integration points. This integration scheme can integrate up to f...
Definition: integral.h:482
3D Gaussian integration class. 4x4x4 integration points. This integration scheme can integrate up to ...
Definition: integral.h:528
bool is_halo() const
Is this element a halo?
Definition: elements.h:1167
Generic class for numerical integration schemes:
Definition: integral.h:49
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
static bool mpi_has_been_initialised()
return true if MPI has been initialised
static OomphCommunicator * communicator_pt()
access to global communicator. This is the oomph-lib equivalent of MPI_COMM_WORLD
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
A general mesh class.
Definition: mesh.h:67
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed, i.e. if we don't have mpi).
Definition: mesh.h:1600
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Vector< GeneralisedElement * > haloed_element_pt(const unsigned &p)
Return vector of haloed elements in this Mesh whose haloing counterpart is held on processor p.
Definition: mesh.h:1779
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Vector< GeneralisedElement * > halo_element_pt(const unsigned &p)
Return vector of halo elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1740
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Telements.h:1130
1D Gaussian integration class for linear "triangular" elements. Two integration points....
Definition: integral.h:640
1D Gaussian integration class for quadratic "triangular" elements. Three integration points....
Definition: integral.h:687
1D Gaussian integration class for cubic "triangular" elements. Four integration points....
Definition: integral.h:734
2D Gaussian integration class for linear triangles. Three integration points. This integration scheme...
Definition: integral.h:820
2D Gaussian integration class for quadratic triangles. Seven integration points. This integration sch...
Definition: integral.h:867
2D Gaussian integration class for cubic triangles. Thirteen integration points. This integration sche...
Definition: integral.h:914
3D Gaussian integration class for tets. Four integration points. This integration scheme can integrat...
Definition: integral.h:1169
3D Gaussian integration class for tets. Eleven integration points. This integration scheme can integr...
Definition: integral.h:1217
3D Gaussian integration class for tets. 45 integration points. This integration scheme can integrate ...
Definition: integral.h:1265
void get_recovered_flux_in_patch(const Vector< ElementWithZ2ErrorEstimator * > &patch_el_pt, const unsigned &num_recovery_terms, const unsigned &num_flux_terms, const unsigned &dim, DenseMatrix< double > *&recovered_flux_coefficient_pt)
Given the vector of elements that make up a patch, the number of recovery and flux terms,...
unsigned nrecovery_terms(const unsigned &dim)
Return number of coefficients for expansion of recovered fluxes for given spatial dimension of elemen...
CombinedErrorEstimateFctPt Combined_error_fct_pt
Function pointer to combined error estimator function.
double get_combined_error_estimate(const Vector< double > &compound_error)
Return a combined error estimate from all compound errors.
void setup_patches(Mesh *&mesh_pt, std::map< Node *, Vector< ElementWithZ2ErrorEstimator * > * > &adjacent_elements_pt, Vector< Node * > &vertex_node_pt)
Setup patches: For each vertex node pointed to by nod_pt, adjacent_elements_pt[nod_pt] contains the p...
unsigned Recovery_order
Order of recovery polynomials.
double Reference_flux_norm
Prescribed reference flux norm.
void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error)
Compute the elemental error measures for a given mesh and store them in a vector.
void shape_rec(const Vector< double > &x, const unsigned &dim, Vector< double > &psi_r)
Recovery shape functions as functions of the global, Eulerian coordinate x of dimension dim....
Integral * integral_rec(const unsigned &dim, const bool &is_q_mesh)
Integation scheme associated with the recovery shape functions must be of sufficiently high order to ...
void doc_flux(Mesh *mesh_pt, const unsigned &num_flux_terms, MapMatrixMixed< Node *, int, double > &rec_flux_map, const Vector< double > &elemental_error, DocInfo &doc_info)
Doc flux and recovered flux.
unsigned & recovery_order()
Access function for order of recovery polynomials.
bool Recovery_order_from_first_element
Bool to indicate if recovery order is to be read out from first element in mesh or set globally.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...