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