projection.h
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$
7 //LIC//
8 //LIC// $LastChangedDate$
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 //Header file for a classes used to represent projectable elements
31 
32 #ifndef OOMPH_PROJECTION_HEADER
33 #define OOMPH_PROJECTION_HEADER
34 
35 
36 #include "mesh.h"
37 #include "problem.h"
38 #include "multi_domain.h"
39 #include "shape.h"
40 #include "element_with_external_element.h"
41 #include "linear_solver.h"
42 
43 // Using CG to solve the projection problem
44 #ifdef OOMPH_HAS_TRILINOS
45 #include "trilinos_solver.h"
46 #endif // #ifdef OOMPH_HAS_TRILINOS
47 #include "iterative_linear_solver.h"
48 
49 // Use a preconditioner for the iterative solver
50 #include "preconditioner.h"
51 #include "general_purpose_preconditioners.h"
52 
53 namespace oomph
54 {
55 
56 //==================================================================
57 ///Template-free Base class for projectable elements
58 //==================================================================
60 {
61 
62  protected:
63 
64  ///\short Enumerated collection to specify which projection problem
65  ///is to be solved.
67 
68  /// Field that is currently being projected
69  unsigned Projected_field;
70 
71  ///Time level we are projecting (0=current values; >0: history values)
73 
74  ///\short When projecting the history values of the nodal coordinates,
75  /// this is the coordinate we're projecting
77 
78  ///\short When projecting the Lagrangain coordinates indicate which
79  ///coordiante is to be projected
81 
82  /// \short Variable to indicate if we're projecting the history values of the
83  /// nodal coordinates (Coordinate) the values themselves (Value), or
84  /// the Lagrangian coordinates in Solid Mechanics problems (Lagrangian)
86 
87  /// \short Bool to know if we do projection or not. If false (the default)
88  /// we solve the element's "real" equations rather than the projection
89  /// equations
91 
92 
93  /// \short Store number of "external" interactions that were assigned to
94  /// the element before doing the projection.
96 
97  /// \short Remember if the element includes external geometric data
98  /// when used in non-projection mode (this is temporarily disabled during the
99  /// projection)
101 
102 
103  /// \short Remember if the element includes external data when used in
104  /// non-projection mode (this is temporarily disabled during the
105  /// projection)
107 
108 
109 public:
110 
111 
112  /// Constructor: Initialise data so that we don't project but solve
113  /// the "normal" equations associated with the element.
114  ProjectableElementBase() : Projected_field(0), Time_level_for_projection(0),
115  Projected_coordinate(0), Projected_lagrangian(0),
116  Projection_type(Value), Do_projection(false),
117  Backup_ninteraction(0), Backup_external_geometric_data(false) {}
118 
119  ///Virtual destructor
121 
122  /// \short Pure virtual function in which the element writer
123  /// must specify the values associated with field fld.
124  /// The information is returned in a vector of pairs which comprise
125  /// the Data object and the value within it, that correspond to field fld.
126  /// E.g. in Taylor Hood elements the fld-th velocities are stored
127  /// at the fld-th value of the nodes; the pressures (the DIM-th
128  /// field) are the DIM-th values at the vertex nodes etc.
129  virtual Vector<std::pair<Data*,unsigned> >
130  data_values_of_field(const unsigned& fld)=0;
131 
132  /// \short Number of fields of the problem, so e.g. for 2D Navier Stokes
133  /// this would be 3 (for the two velocities and one pressure)
134  virtual unsigned nfields_for_projection()=0;
135 
136  /// \short Number of history values to be stored for fld-th field
137  /// (includes current value!)
138  virtual unsigned nhistory_values_for_projection(const unsigned &fld)=0;
139 
140  ///\short Number of history values to be stored when projecting
141  /// the history values of the nodal coordinates (includes current value!)
142  virtual unsigned nhistory_values_for_coordinate_projection()=0;
143 
144  /// \short Return number of values (pinned or not) associated with
145  /// field fld within the element. This must correspond to the
146  /// number of shape functions returned in jacobian_and_shape_of_field(...).
147  virtual unsigned nvalue_of_field(const unsigned &fld)=0;
148 
149  /// \short Return local equation numbers associated with value ivalue
150  /// of field fld within the element.
151  virtual int local_equation(const unsigned &fld,const unsigned &ivalue)=0;
152 
153  /// \short Return Jacobian of mapping and the shape functions associated
154  /// with field fld. The number of shape functions must match the
155  /// number of values specified in nvalue_of_field(...). For
156  /// Lagrange-type interpolations the shape functinos are simply
157  /// the "normal" nodal shape functions; if the element contains
158  /// internal Data that is not associated with shape functions,
159  /// simply set the corresonding shape function to 1.
160  virtual double jacobian_and_shape_of_field
161  (const unsigned &fld, const Vector<double> &s, Shape &psi)=0;
162 
163  /// \short Return the fld-th field at local coordinates s
164  /// at time-level time (time=0: current value; time>0: history values)
165  virtual double get_field
166  (const unsigned &time, const unsigned &fld,const Vector<double> &s)=0;
167 
168 };
169 
170 
171 
172 
173 
174 //=====================================================================
175 /// Wrapper class for projectable elements. Adds "projectability"
176 /// to the underlying ELEMENT.
177 //=====================================================================
178 template<class ELEMENT>
179 class ProjectableElement : public virtual ELEMENT,
180  public virtual ProjectableElementBase,
181  public virtual ElementWithExternalElement
182 {
183 
184  protected:
185 
186  /// Overloaded version of fill_in_contribution_to_residuals
187  void fill_in_contribution_to_residuals(Vector<double> &residuals)
188  {
189  //Do projection
190  if (Do_projection)
191  {
192  this->residual_for_projection
193  (residuals,GeneralisedElement::Dummy_matrix, 0);
194  }
195  //solve problem normally
196  else
197  {
198  ELEMENT::fill_in_contribution_to_residuals(residuals);
199  }
200  }
201 
202  /// \short Function to describe the local dofs of the element. The ostream
203  /// specifies the output stream to which the description
204  /// is written; the string stores the currently
205  /// assembled output that is ultimately written to the
206  /// output stream by Data::describe_dofs(...); it is typically
207  /// built up incrementally as we descend through the
208  /// call hierarchy of this function when called from
209  /// Problem::describe_dofs(...)
210  void describe_local_dofs(std::ostream& out,
211  const std::string& current_string) const
212  {
213  ElementWithExternalElement::describe_local_dofs(out,current_string);
214  ELEMENT::describe_local_dofs(out,current_string);
215  }
216 
217  ///Overloaded version of fill_in_contribution_to_jacobian
218  void fill_in_contribution_to_jacobian(Vector<double> &residuals,
219  DenseMatrix<double> &jacobian)
220  {
221  //Do projection
222  if (Do_projection)
223  {
224  this->residual_for_projection(residuals,jacobian,1);
225  }
226  else
227  {
228  ELEMENT::fill_in_contribution_to_jacobian(residuals,jacobian);
229  }
230  }
231 
232 
233  public:
234 
235 
236  /// \short Constructor [this was only required explicitly
237  /// from gcc 4.5.2 onwards...]
239 
240  /// \short Residual for the projection step. Flag indicates if we
241  /// want the Jacobian (1) or not (0). Virtual so it can be
242  /// overloaded if necessary
243  virtual void residual_for_projection(Vector<double> &residuals,
244  DenseMatrix<double> &jacobian,
245  const unsigned& flag)
246  {
247  //Am I a solid element
248  SolidFiniteElement* solid_el_pt =
249  dynamic_cast<SolidFiniteElement*>(this);
250 
251  unsigned n_dim=dim();
252 
253  //Allocate storage for local coordinates
254  Vector<double> s(n_dim);
255 
256  //Current field
257  unsigned fld=Projected_field;
258 
259  //Number of nodes
260  const unsigned n_node = this->nnode();
261  //Number of positional dofs
262  const unsigned n_position_type = this->nnodal_position_type();
263 
264  //Number of dof for current field
265  const unsigned n_value=nvalue_of_field(fld);
266 
267  //Set the value of n_intpt
268  const unsigned n_intpt = integral_pt()->nweight();
269 
270  //Loop over the integration points
271  for(unsigned ipt=0;ipt<n_intpt;ipt++)
272  {
273  // Get the local coordinates of Gauss point
274  for(unsigned i=0;i<n_dim;i++) s[i] = integral_pt()->knot(ipt,i);
275 
276  //Get the integral weight
277  double w = integral_pt()->weight(ipt);
278 
279  // Find same point in base mesh using external storage
280  FiniteElement* other_el_pt=0;
281  other_el_pt=this->external_element_pt(0,ipt);
282  Vector<double> other_s(n_dim);
283  other_s=this->external_element_local_coord(0,ipt);
284 
285  ProjectableElement<ELEMENT>* cast_el_pt =
286  dynamic_cast<ProjectableElement<ELEMENT>*>(other_el_pt);
287 
288  //Storage for the local equation and local unknown
289  int local_eqn=0, local_unknown=0;
290 
291  //Now set up the three different projection problems
292  switch(Projection_type)
293  {
294  case Lagrangian:
295  {
296  //If we don't have a solid element there is a problem
297  if(solid_el_pt==0)
298  {
299  throw OomphLibError(
300  "Trying to project Lagrangian coordinates in non-SolidElement\n",
301  OOMPH_CURRENT_FUNCTION,
302  OOMPH_EXCEPTION_LOCATION);
303  }
304 
305  //Find the position shape functions
306  Shape psi(n_node,n_position_type);
307  //Get the position shape functions
308  this->shape(s,psi);
309  //Get the jacobian
310  double J = this->J_eulerian(s);
311 
312  //Premultiply the weights and the Jacobian
313  double W = w*J;
314 
315  //Get the value of the current position of the 0th coordinate
316  //in the current element
317  //at the current time level (which is the unkown)
318  double interpolated_xi_proj = this->interpolated_x(s,0);
319 
320  //Get the Lagrangian position in the other element
321  double interpolated_xi_bar=
322  dynamic_cast<SolidFiniteElement*>(cast_el_pt)
323  ->interpolated_xi(other_s,Projected_lagrangian);
324 
325  //Now loop over the nodes and position dofs
326  for(unsigned l=0;l<n_node;++l)
327  {
328  //Loop over position unknowns
329  for(unsigned k=0;k<n_position_type;++k)
330  {
331  //The local equation is going to be the positional local equation
332  local_eqn = solid_el_pt->position_local_eqn(l,k,0);
333 
334  //Now assemble residuals
335  if(local_eqn >= 0)
336  {
337  //calculate residuals
338  residuals[local_eqn] +=
339  (interpolated_xi_proj - interpolated_xi_bar)*psi(l,k)*W;
340 
341  //Calculate the jacobian
342  if(flag==1)
343  {
344  for(unsigned l2=0;l2<n_node;l2++)
345  {
346  //Loop over position dofs
347  for(unsigned k2=0;k2<n_position_type;k2++)
348  {
349  local_unknown =
350  solid_el_pt->position_local_eqn(l2,k2,0);
351  if(local_unknown >= 0)
352  {
353  jacobian(local_eqn,local_unknown)
354  += psi(l2,k2)*psi(l,k)*W;
355  }
356  }
357  }
358  } //end of jacobian
359  }
360  }
361  }
362  } //End of Lagrangian coordinate case
363 
364  break;
365 
366  //Now the coordinate history case
367  case Coordinate:
368  {
369  //Find the position shape functions
370  Shape psi(n_node,n_position_type);
371  //Get the position shape functions
372  this->shape(s,psi);
373  //Get the jacobian
374  double J = this->J_eulerian(s);
375 
376  //Premultiply the weights and the Jacobian
377  double W = w*J;
378 
379  //Get the value of the current position in the current element
380  //at the current time level (which is the unkown)
381  double interpolated_x_proj = 0.0;
382  //If we are a solid element read it out directly from the data
383  if(solid_el_pt!=0)
384  {
385  interpolated_x_proj = this->interpolated_x(s,Projected_coordinate);
386  }
387  //Otherwise it's the field value at the current time
388  else
389  {
390  interpolated_x_proj = this->get_field(0,fld,s);
391  }
392 
393  //Get the position in the other element
394  double interpolated_x_bar=
395  cast_el_pt->interpolated_x(Time_level_for_projection,
396  other_s,Projected_coordinate);
397 
398  //Now loop over the nodes and position dofs
399  for(unsigned l=0;l<n_node;++l)
400  {
401  //Loop over position unknowns
402  for(unsigned k=0;k<n_position_type;++k)
403  {
404  //If I'm a solid element
405  if(solid_el_pt!=0)
406  {
407  //The local equation is going to be the positional local equation
408  local_eqn =
409  solid_el_pt->position_local_eqn(l,k,Projected_coordinate);
410  }
411  //Otherwise just pick the local unknown of a field to
412  //project into
413  else
414  {
415  //Complain if using generalised position types
416  //but this is not a solid element, because the storage
417  //is then not clear
418  if(n_position_type!=1)
419  {
420  throw OomphLibError(
421  "Trying to project generalised positions not in SolidElement\n",
422  OOMPH_CURRENT_FUNCTION,
423  OOMPH_EXCEPTION_LOCATION);
424  }
425  local_eqn = local_equation(fld,l);
426  }
427 
428  //Now assemble residuals
429  if(local_eqn >= 0)
430  {
431  //calculate residuals
432  residuals[local_eqn] +=
433  (interpolated_x_proj - interpolated_x_bar)*psi(l,k)*W;
434 
435  //Calculate the jacobian
436  if(flag==1)
437  {
438  for(unsigned l2=0;l2<n_node;l2++)
439  {
440  //Loop over position dofs
441  for(unsigned k2=0;k2<n_position_type;k2++)
442  {
443  //If I'm a solid element
444  if(solid_el_pt!=0)
445  {
446  local_unknown = solid_el_pt
447  ->position_local_eqn(l2,k2,Projected_coordinate);
448  }
449  else
450  {
451  local_unknown = local_equation(fld,l2);
452  }
453 
454  if(local_unknown >= 0)
455  {
456  jacobian(local_eqn,local_unknown)
457  += psi(l2,k2)*psi(l,k)*W;
458  }
459  }
460  }
461  } //end of jacobian
462  }
463  }
464  }
465  } //End of coordinate case
466  break;
467 
468  //Now the value case
469  case Value:
470  {
471  //Field shape function
472  Shape psi(n_value);
473 
474  //Calculate jacobian and shape functions for that field
475  double J=jacobian_and_shape_of_field(fld,s,psi);
476 
477  //Premultiply the weights and the Jacobian
478  double W = w*J;
479 
480  //Value of field in current element at current time level
481  //(the unknown)
482  unsigned now=0;
483  double interpolated_value_proj = this->get_field(now,fld,s);
484 
485  //Value of the interpolation of element located in base mesh
486  double interpolated_value_bar =
487  cast_el_pt->get_field(Time_level_for_projection,fld,other_s);
488 
489  //Loop over dofs of field fld
490  for(unsigned l=0;l<n_value;l++)
491  {
492  local_eqn = local_equation(fld,l);
493  if(local_eqn >= 0)
494  {
495  //calculate residuals
496  residuals[local_eqn] +=
497  (interpolated_value_proj - interpolated_value_bar)*psi[l]*W;
498 
499  //Calculate the jacobian
500  if(flag==1)
501  {
502  for(unsigned l2=0;l2<n_value;l2++)
503  {
504  local_unknown = local_equation(fld,l2);
505  if(local_unknown >= 0)
506  {
507  jacobian(local_eqn,local_unknown)
508  += psi[l2]*psi[l]*W;
509  }
510  }
511  } //end of jacobian
512  }
513  }
514  }
515  break;
516 
517  default:
518  throw OomphLibError(
519  "Wrong type specificied in Projection_type. This should never happen\n",
520  OOMPH_CURRENT_FUNCTION,
521  OOMPH_EXCEPTION_LOCATION);
522  } //End of the switch statement
523 
524  }//End of loop over ipt
525 
526  }//End of residual_for_projection function
527 
528 
529  /// \short Use Eulerian coordinates for matching in locate_zeta
530  /// when doing projection
531  double zeta_nodal(const unsigned &n, const unsigned &k,
532  const unsigned &i) const
533  {
534  if (Do_projection)
535  {
536  return nodal_position_gen(n,k,i);
537  }
538  else
539  {
540  return ELEMENT::zeta_nodal(n,k,i);
541  }
542  }
543 
544 
545 
546  /// \short Backup the element's state and
547  /// switch it to projection mode.
549  {
550  //Backup number of interaction
551  Backup_ninteraction= ninteraction();
552 
553  //Backup flag for inclusion of geometric data
554  if (add_external_geometric_data())
555  {
557  }
558  else
559  {
561  }
562 
563  //Backup flag for inclusion of interaction data
564  if (add_external_interaction_data())
565  {
567  }
568  else
569  {
571  }
572 
573  //Actions to enable projection
574  Do_projection=true;
575  ignore_external_geometric_data();
576  ignore_external_interaction_data();
577  set_ninteraction(1);
578  }
579 
580  ///\short Helper function to restore the element to the state
581  /// it was in before we entered the projection mode and switch off
582  /// projection mode.
584  {
585  //Restore number of interaction
586  set_ninteraction(Backup_ninteraction);
587 
588  //Restore geometric data
590  {
591  include_external_geometric_data();
592  }
593  else
594  {
595  ignore_external_geometric_data();
596  }
597 
598  //Restore interaction data
600  {
601  include_external_interaction_data();
602  }
603  else
604  {
605  ignore_external_interaction_data();
606  }
607 
608  Do_projection=false;
609  }
610 
611 
612  ///Project (history values of) coordintes
614 
615  /// Project (history values of) values
617 
618  /// Project (current and only values of) lagrangian coordinates
620 
621 
622  ///Field that is currently being projected
623  unsigned& projected_field()
624  {return Projected_field;}
625 
626 ///Which history value are we projecting?
628  {return Time_level_for_projection;}
629 
630  /// \short When projecting the history values of the nodal coordinates,
631  /// this is the coordinate we're projecting
633  {return Projected_coordinate;}
634 
635  /// \short When projecting the Lagrangian coordinates this is
636  /// the coordinate that is being projected
638  {return Projected_lagrangian;}
639 
640 
641 };//End of class
642 
643 
644 
645 
646 
647 //=======================================================================
648 /// Face geometry for element is the same as that for the underlying
649 /// wrapped element
650 //=======================================================================
651 template<class ELEMENT>
652 class FaceGeometry<ProjectableElement<ELEMENT> >
653 : public virtual FaceGeometry<ELEMENT>
654  {
655  public:
656  FaceGeometry() : FaceGeometry<ELEMENT>() {}
657  };
658 
659 // Forward definition of the friends of the class
660 
661 // The RefineableTriangleMesh
662 //template<class FRIEND_PROJECTABLE_ELEMENT>
663 //class RefineableTriangleMesh;
664 
665 // The RefineableTetgenMesh
666 //template<class FRIEND_PROJECTABLE_ELEMENT>
667 //class RefineableTetgenMesh;
668 
669 // The BackupMeshForProjection
670 //template<class FRIEND_PROJECTABLE_ELEMENT>
671 //class BackupMeshForProjection;
672 
673 //=======================================================================
674 /// Projection problem. This is created during the adaptation
675 /// of unstructured meshes and it is assumed that no boundary conditions
676 /// have been set. If they have, they will be unset during the projection
677 /// and must be reset afterwards.
678 //=======================================================================
679 template<class PROJECTABLE_ELEMENT>
680 class ProjectionProblem : public virtual Problem
681 {
682 
683  // The classes are friend whether the templated element of the
684  // friend class is the same or not as the templated element of the
685  // ProjectionProblem class
686  template<class FRIEND_PROJECTABLE_ELEMENT> friend class RefineableTriangleMesh;
687  template<class FRIEND_PROJECTABLE_ELEMENT> friend class RefineableTetgenMesh;
688  template<class FRIEND_PROJECTABLE_ELEMENT> friend class BackupMeshForProjection;
689  template<class FRIEND_PROJECTABLE_ELEMENT> friend class RefineableGmshTetMesh;
690 
691  // The classes are friend only when the templated element of the
692  // friend class matches the templated element of the
693  // ProjectionProblem class
694  // friend class RefineableTriangleMesh<class PROJECTABLE_ELEMENT>;
695  // friend class RefineableTetgenMesh<class PROJECTABLE_ELEMENT>;
696  // friend class BackupMeshForProjection<class PROJECTABLE_ELEMENT>;
697 
698  public:
699 
700  /// Suppress all output during projection phases
702  {
703  Output_during_projection_suppressed=true;
704  }
705 
706  /// Undo suppression of all output during projection phases
708  {
709  Output_during_projection_suppressed=false;
710  }
711 
712  /// Return the value of the flag about using an iterative solver for
713  /// projection
715  {return Use_iterative_solver_for_projection;}
716 
717  /// Enables the use of an iterative solver for projection
719  {Use_iterative_solver_for_projection=true;}
720 
721  /// Disbales the use of an iterative solver for projection
723  {Use_iterative_solver_for_projection=false;}
724 
725  ///\short Project from base into the problem's own mesh.
726  void project(Mesh* base_mesh_pt, const bool& dont_project_positions=false)
727  {
728  // Use an iterative solver?
729  if (Use_iterative_solver_for_projection)
730  {
731  // If oomph-lib has Trilinos installed then use the CG version
732  // from Trilinos, otherwise use oomph-lib's own CG implementation
733 #ifdef OOMPH_HAS_TRILINOS
734  // Check whether the problem is distributed?
735  if (MPI_Helpers::mpi_has_been_initialised())
736  {
737  // Create a Trilinos Solver
738  Iterative_solver_projection_pt = new TrilinosAztecOOSolver;
739  // Select CG as the linear solver
740  dynamic_cast<TrilinosAztecOOSolver*>(Iterative_solver_projection_pt)->
741  solver_type() = TrilinosAztecOOSolver::CG;
742  }
743  else
744  {
745  // Use CG to solve the projection problem
746  Iterative_solver_projection_pt = new CG<CRDoubleMatrix>;
747  }
748 
749  // Create the preconditioner
750  Preconditioner_projection_pt = new MatrixBasedDiagPreconditioner();
751  // Set the preconditioner for the solver
752  Iterative_solver_projection_pt->preconditioner_pt() =
753  Preconditioner_projection_pt;
754 
755  // Set CG as the linear solver
756  Problem::linear_solver_pt() = Iterative_solver_projection_pt;
757 
758 #else
759  // Check whether the problem is distributed?
760  if (!(MPI_Helpers::mpi_has_been_initialised()))
761  {
762  // If we did not installed Trilinos and the problem is not
763  // distributed then we can use a (serial) preconditioned
764  // iterative solver, otherwise, if we did not installed Trilinos
765  // but the problem is distributed then we cannot use a
766  // preconditioned iterative solver. Matrix multiplication in a
767  // distributed environment is only performed by Trilinos. We
768  // then use a direct solver for the projection problem.
769 
770  // Use CG to solve the projection problem
771  Iterative_solver_projection_pt = new CG<CRDoubleMatrix>;
772 
773  // Create the preconditioner
774  Preconditioner_projection_pt = new MatrixBasedDiagPreconditioner();
775  // Set the preconditioner for the solver
776  Iterative_solver_projection_pt->preconditioner_pt() =
777  Preconditioner_projection_pt;
778 
779  // Set CG as the linear solver
780  Problem::linear_solver_pt() = Iterative_solver_projection_pt;
781  }
782  else
783  {
784  // Use a direct solver. Do nothing
785  }
786 
787 #endif
788 
789  } // if (Use_iterative_solver_for_projection)
790 
791  // Backup verbosity in Newton solve status
792  bool shut_up_in_newton_solve_backup=Shut_up_in_newton_solve;
793 
794  // Disable documentation of solve times
795  bool backed_up_doc_time_enabled=linear_solver_pt()->is_doc_time_enabled();
796  if (Output_during_projection_suppressed)
797  {
798  linear_solver_pt()->disable_doc_time();
799  }
800 
801  //Display stats
802  unsigned n_element = Problem::mesh_pt()->nelement();
803  unsigned n_element1=base_mesh_pt->nelement();
804  unsigned n_node = Problem::mesh_pt()->nnode();
805  unsigned n_node1=base_mesh_pt->nnode();
806  if (!Output_during_projection_suppressed)
807  {
808  oomph_info <<"\n=============================\n";
809  oomph_info << "Base mesh has " << n_element1 << " elements\n";
810  oomph_info << "Target mesh has " << n_element << " elements\n";
811  oomph_info << "Base mesh has " << n_node1 << " nodes\n";
812  oomph_info << "Target mesh has " << n_node << " nodes\n";
813  oomph_info <<"=============================\n\n";
814 
815  }
816  else
817  {
818  // Make Newton solver shut up too
819  disable_info_in_newton_solve();
820  }
821 
822 
823  if (n_element==0)
824  {
825  oomph_info
826  << "Very odd -- no elements in target mesh; "
827  << " not doing anything in ProjectionProblem::project()\n";
828  return;
829  }
830 
831 #ifdef PARANOID
832  unsigned nnod=Problem::mesh_pt()->nnode();
833  if (nnod==0)
834  {
835  std::ostringstream error_stream;
836  error_stream
837  << "Mesh has no nodes! Please populate the Node_pt vector\n"
838  << "otherwise projection won't work!\n";
839  throw OomphLibError(error_stream.str(),
840  OOMPH_CURRENT_FUNCTION,
841  OOMPH_EXCEPTION_LOCATION);
842  }
843 #endif
844 
845  // How many fields do we have to project?
846  unsigned n_fields = dynamic_cast<PROJECTABLE_ELEMENT*>
847  (Problem::mesh_pt()->element_pt(0))->nfields_for_projection();
848 
849  // Spatial dimension of the problem
850  unsigned n_dim = Problem::mesh_pt()->node_pt(0)->ndim();
851 
852  // Default number of history values
853  unsigned n_history_values=0;
854 
855  // Set things up for coordinate projection
856  for(unsigned e=0;e<n_element;e++)
857  {
858  PROJECTABLE_ELEMENT * el_pt =
859  dynamic_cast<PROJECTABLE_ELEMENT*>
860  (Problem::mesh_pt()->element_pt(e));
861 
862  // Switch to projection
863  el_pt->enable_projection();
864  }
865 
866  // Switch elements in base mesh to projection mode (required
867  // to switch to use of Eulerian coordinates when identifying
868  // corresponding points in the two meshes)
869  for(unsigned e=0;e<n_element1;e++)
870  {
871  PROJECTABLE_ELEMENT * el_pt =
872  dynamic_cast<PROJECTABLE_ELEMENT*>
873  (base_mesh_pt->element_pt(e));
874 
875  // Switch to projection
876  el_pt->enable_projection();
877  }
878 
879 
880  // Set up multi domain interactions so we can locate the
881  // values in the base mesh.
882  // Note that it's important to switch elements to projection
883  // mode first to ensure that matching is done based on Eulerian
884  // rather than Lagrangian coordinates if pseudo-solid elements
885  // are used.
886  double t_start=TimingHelpers::timer();
887  Multi_domain_functions::setup_multi_domain_interaction
888  <PROJECTABLE_ELEMENT>(this,Problem::mesh_pt(),base_mesh_pt);
889  if (!Output_during_projection_suppressed)
890  {
891  oomph_info <<"CPU for setup of multi-domain interaction for projection: "
892  << TimingHelpers::timer()-t_start << std::endl;
893  }
894  t_start=TimingHelpers::timer();
895 
896 
897  //Let us first pin every degree of freedom
898  //We shall unpin selected dofs for each different projection problem
899  this->pin_all();
900 
901  if (!dont_project_positions)
902  {
903 
904  //------------------Project coordinates first------------------------
905  //If we have a solid element then we should also project Lagrangian
906  //coordinates, but we can use the storage that MUST be provided for
907  //the unknown positions for this.
908  //If we can cast the first element of the mesh to a solid element,
909  //then assume that we have solid elements
910  if(dynamic_cast<SolidFiniteElement*>(
911  Problem::mesh_pt()->element_pt(0)))
912  {
913  // Store current positions
914  this->store_positions();
915 
916  //There are no history values for the Lagrangian coordinates
917  //Set coordinate 0 for projection
918  this->set_coordinate_for_projection(0);
919  this->unpin_dofs_of_coordinate(0);
920 
921  //Loop over the Lagrangian coordinate
922  const unsigned n_lagrangian =
923  dynamic_cast<SolidNode*>(Problem::mesh_pt()->node_pt(0))->nlagrangian();
924 
925  //Now loop over the lagrangian coordinates
926  for(unsigned i=0;i<n_lagrangian;++i)
927  {
928 
929  if (!Output_during_projection_suppressed)
930  {
931  oomph_info <<"\n\n=============================================\n";
932  oomph_info << "Projecting values for Lagrangian coordinate " << i
933  << std::endl;
934  oomph_info <<"=============================================\n\n";
935  }
936 
937  //Set the coordinate for projection
938  this->set_lagrangian_coordinate_for_projection(i);
939 
940  //Assign equation number
941  unsigned ndof_tmp=assign_eqn_numbers();
942  if (!Output_during_projection_suppressed)
943  {
944  oomph_info <<
945  "Number of equations for projection of Lagrangian coordinate "
946  << " : "<< ndof_tmp <<std::endl << std::endl;
947  }
948 
949 
950  if (Problem_is_nonlinear)
951  {
952  std::ostringstream error_stream;
953  error_stream
954  << "Solid mechanics problems will break if Problem_is_nonlinear is\n"
955  << "set to true in the projection problem because we're using the\n "
956  << "actual nodal positions to store the values of the Lagrangian\n"
957  << "coords. There shouldn't be any need for \n"
958  << "Problem_is_nonlinear=true anyway, apart from debugging in \n"
959  << "which case you now know why this case will break!\n";
960  OomphLibWarning(error_stream.str(),
961  OOMPH_CURRENT_FUNCTION,
962  OOMPH_EXCEPTION_LOCATION);
963  }
964 
965 
966  //Projection and interpolation
967  Problem::newton_solve();
968 
969  //Move values back into Lagrangian coordinate for all nodes
970  unsigned n_node = Problem::mesh_pt()->nnode();
971  for(unsigned n=0;n<n_node;++n)
972  {
973  //Cast it to a solid node
974  SolidNode* solid_node_pt =
975  dynamic_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
976  //Now find number of types
977  const unsigned n_position_type = solid_node_pt->nposition_type();
978  //Find number of lagrangian types
979  const unsigned n_lagrangian_type = solid_node_pt->nlagrangian_type();
980 
981  //If these are not the same, throw an error
982  if(n_position_type != n_lagrangian_type)
983  {
984  std::ostringstream error_stream;
985  error_stream
986  << "The number of generalised position dofs "
987  << n_position_type
988  << "\n not the same as the number of generalised lagrangian dofs "
989  << n_lagrangian_type << ".\n"
990  << "Projection cannot be done at the moment, sorry.\n";
991 
992  throw OomphLibError(error_stream.str(),
993  OOMPH_CURRENT_FUNCTION,
994  OOMPH_EXCEPTION_LOCATION);
995  }
996 
997  //Now transfer the information across
998  //from the first coordinate which was used during the projection
999  for(unsigned k=0;k<n_position_type;++k)
1000  {
1001  solid_node_pt->xi_gen(k,i) = solid_node_pt->x_gen(k,0);
1002  //Reset real values so that the Jacobians are correctly
1003  //computed next time around
1004  solid_node_pt->x_gen(k,0) = Solid_backup[n](k,0);
1005  }
1006  }
1007  } //End of loop over lagrangian coordinates
1008 
1009  //Now repin the dofs
1010  this->pin_dofs_of_coordinate(0);
1011 
1012  //Now project the position histories
1013 
1014  //Check number of history values for coordinates
1015  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1016  (Problem::mesh_pt()->element_pt(0))->
1018 
1019  //Projection the coordinates only if there are history values
1020  if(n_history_values > 1)
1021  {
1022  for (unsigned i=0;i<n_dim;i++)
1023  {
1024  if (!Output_during_projection_suppressed)
1025  {
1026  oomph_info <<"\n\n=============================================\n";
1027  oomph_info << "Projecting history values for coordinate " << i
1028  << std::endl;
1029  oomph_info <<"=============================================\n\n";
1030  }
1031 
1032  //Set the coordinate for projection
1033  this->set_coordinate_for_projection(i);
1034  this->unpin_dofs_of_coordinate(i);
1035 
1036  // Loop over number of history values, beginning with the latest one.
1037  // Don't deal with current time.
1038  for (unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1039  {
1040  unsigned time_level=h_tim-1;
1041 
1042  //Set time_level we are dealing with
1043  this->set_time_level_for_projection(time_level);
1044 
1045  //Assign equation number
1046  unsigned ndof_tmp=assign_eqn_numbers();
1047  if (!Output_during_projection_suppressed)
1048  {
1049  oomph_info << "Number of equations for projection of coordinate "
1050  << i << " at time level " << time_level
1051  << " : "<< ndof_tmp <<std::endl << std::endl;
1052  }
1053 
1054  //Projection and interpolation
1055  Problem::newton_solve();
1056 
1057  //Move values back into history value of coordinate
1058  unsigned n_node = Problem::mesh_pt()->nnode();
1059  for(unsigned n=0;n<n_node;++n)
1060  {
1061  //Cache the pointer to the node
1062  Node* nod_pt = Problem::mesh_pt()->node_pt(n);
1063  //Find the number of generalised dofs
1064  const unsigned n_position_type = nod_pt->nposition_type();
1065  //Now copy all back
1066  for(unsigned k=0;k<n_position_type;++k)
1067  {
1068  nod_pt->x_gen(time_level,k,i) = nod_pt->x_gen(0,k,i);
1069  //Reset real values so that the Jacobians are computed
1070  //correctly next time around
1071  nod_pt->x_gen(0,k,i) = Solid_backup[n](k,i);
1072  }
1073  }
1074  }
1075  //Repin the positions
1076  this->pin_dofs_of_coordinate(i);
1077  }
1078  } //End of history value projection
1079  } //End of SolidElement case
1080 
1081  //Now for non solid elements, we are going to hijack the
1082  //first value as storage to be used for the projection of the history
1083  //values
1084  else
1085  {
1086  // Prepare for projection in value 0
1087  this->set_current_field_for_projection(0);
1088  this->unpin_dofs_of_field(0);
1089 
1090 #ifdef PARANOID
1091 
1092  // The machinery used below assumes that field 0 can actually
1093  // hold the coordinates i.e. that the field is interpolated
1094  // isoparametrically! The method will fail if there are no values
1095  // stored at the nodes. Currently there are no examples of that -- the
1096  // problem would only arise for elements that have all their fields
1097  // represented by internal data. Will fix this if/when
1098  // such a case arises...
1099  unsigned n_element = Problem::mesh_pt()->nelement();
1100  for(unsigned e=0;e<n_element;e++)
1101  {
1102  FiniteElement* el_pt=Problem::mesh_pt()->finite_element_pt(e);
1103  unsigned nnod=el_pt->nnode();
1104  for (unsigned j=0;j<nnod;j++)
1105  {
1106  Node* nod_pt=el_pt->node_pt(j);
1107  if (nod_pt->nvalue()==0)
1108  {
1109  std::ostringstream error_stream;
1110  error_stream << "Node at ";
1111  unsigned n=nod_pt->ndim();
1112  for (unsigned i=0;i<n;i++)
1113  {
1114  error_stream << nod_pt->x(i) << " ";
1115  }
1116  error_stream
1117  << "\nhas no values. Projection (of coordinates) doesn't work\n"
1118  << "for such cases (at the moment), sorry! Send us the code\n"
1119  << "where the problem arises and we'll try to implement this\n"
1120  << "(up to now unnecessary) capability...\n";
1121  throw OomphLibError(error_stream.str(),
1122  OOMPH_CURRENT_FUNCTION,
1123  OOMPH_EXCEPTION_LOCATION);
1124  }
1125  }
1126  }
1127 
1128 #endif
1129 
1130  //Check number of history values for coordinates
1131  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1132  (Problem::mesh_pt()->element_pt(0))->
1134 
1135  //Projection the coordinates only if there are history values
1136  if(n_history_values > 1)
1137  {
1138  for (unsigned i=0;i<n_dim;i++)
1139  {
1140  if (!Output_during_projection_suppressed)
1141  {
1142  oomph_info <<"\n\n=============================================\n";
1143  oomph_info << "Projecting history values for coordinate " << i
1144  << std::endl;
1145  oomph_info <<"=============================================\n\n";
1146  }
1147 
1148  //Set the coordinate for projection
1149  this->set_coordinate_for_projection(i);
1150 
1151  // Loop over number of history values, beginning with the latest one.
1152  // Don't deal with current time.
1153  for (unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1154  {
1155  unsigned time_level=h_tim-1;
1156 
1157  //Set time_level we are dealing with
1158  this->set_time_level_for_projection(time_level);
1159 
1160  //Assign equation number
1161  unsigned ndof_tmp=assign_eqn_numbers();
1162  if (!Output_during_projection_suppressed)
1163  {
1164  oomph_info << "Number of equations for projection of coordinate "
1165  << i << " at time level " << time_level
1166  << " : "<< ndof_tmp <<std::endl << std::endl;
1167  }
1168 
1169  //Projection and interpolation
1170  Problem::newton_solve();
1171 
1172  //Move values back into history value of coordinate
1173  unsigned n_element = Problem::mesh_pt()->nelement();
1174  for(unsigned e=0;e<n_element;e++)
1175  {
1176  PROJECTABLE_ELEMENT * new_el_pt =
1177  dynamic_cast<PROJECTABLE_ELEMENT*>
1178  (Problem::mesh_pt()->element_pt(e));
1179 
1180  Vector<std::pair<Data*,unsigned> >
1181  data=new_el_pt->data_values_of_field(0);
1182 
1183  unsigned d_size=data.size();
1184  for(unsigned d=0;d<d_size;d++)
1185  {
1186  //Replace as coordinates
1187  double coord=data[d].first->value(0,0);
1188  dynamic_cast<Node*>(data[d].first)->x(time_level,i) = coord;
1189  }
1190  }
1191  }
1192  }
1193  } //End of history value projection
1194 
1195  //Repin the dofs for field 0
1196  this->pin_dofs_of_field(0);
1197 
1198  } //End of non-SolidElement case
1199 
1200 
1201  } // end if for projection of coordinates
1202 
1203  //Disable projection of coordinates
1204  for(unsigned e=0;e<n_element;e++)
1205  {
1206  PROJECTABLE_ELEMENT * el_pt =
1207  dynamic_cast<PROJECTABLE_ELEMENT*>
1208  (Problem::mesh_pt()->element_pt(e));
1209 
1210  el_pt->set_project_values();
1211  }
1212 
1213  //Loop over fields
1214  for (unsigned fld=0; fld<n_fields ;fld++)
1215  {
1216  //Let us first pin every degree of freedom
1217  //We shall unpin selected dofs for each different projection problem
1218  this->pin_all();
1219 
1220  //Do actions for this field
1221  this->set_current_field_for_projection(fld);
1222  this->unpin_dofs_of_field(fld);
1223 
1224  //Check number of history values
1225  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1226  (Problem::mesh_pt()->element_pt(0))->nhistory_values_for_projection(fld);
1227 
1228  //Loop over number of history values
1229  //Beginning with the latest one
1230  for (unsigned h_tim=n_history_values;h_tim>0;h_tim--)
1231  {
1232  unsigned time_level=h_tim-1;
1233  if (!Output_during_projection_suppressed)
1234  {
1235  oomph_info <<"\n=========================================\n";
1236  oomph_info << "Projecting field " << fld << " at time level "
1237  << time_level<<std::endl;
1238  oomph_info << "========================================\n";
1239  }
1240 
1241  //Set time_level we are dealing with
1242  this->set_time_level_for_projection(time_level);
1243 
1244  //Assign equation number
1245  unsigned ndof_tmp=assign_eqn_numbers();
1246  if (!Output_during_projection_suppressed)
1247  {
1248  oomph_info << "Number of equations for projection of field "
1249  << fld << " at time level " << time_level
1250  << " : "<< ndof_tmp <<std::endl << std::endl;
1251  }
1252 
1253  //Projection and interpolation
1254  Problem::newton_solve();
1255 
1256  // Move computed values into the required time-level (not needed
1257  // for current values which are done last -- they simply
1258  // stay where they are)
1259  if(time_level!=0)
1260  {
1261  for(unsigned e=0;e<n_element;e++)
1262  {
1263  PROJECTABLE_ELEMENT * new_el_pt =
1264  dynamic_cast<PROJECTABLE_ELEMENT*>
1265  (Problem::mesh_pt()->element_pt(e));
1266 
1267  Vector<std::pair<Data*,unsigned> >
1268  data=new_el_pt->data_values_of_field(fld);
1269 
1270  unsigned d_size=data.size();
1271  for(unsigned d=0;d<d_size;d++)
1272  {
1273  //Move into time level
1274  double c_value=data[d].first->value(0,data[d].second);
1275  data[d].first->set_value(time_level,data[d].second,c_value);
1276  }
1277  }
1278  }
1279  } //End of loop over time levels
1280 
1281  } //End of loop over fields
1282 
1283 
1284  //Reset parameters of external storage and interactions
1285  for(unsigned e=0;e<n_element;e++)
1286  {
1287  PROJECTABLE_ELEMENT * new_el_pt =
1288  dynamic_cast<PROJECTABLE_ELEMENT*>
1289  (Problem::mesh_pt()->element_pt(e));
1290 
1291  // Flush information associated with the external elements
1292  new_el_pt->flush_all_external_element_storage();
1293 
1294  new_el_pt->disable_projection();
1295  }
1296 
1297  for(unsigned e=0;e<n_element1;e++)
1298  {
1299  PROJECTABLE_ELEMENT * el_pt =
1300  dynamic_cast<PROJECTABLE_ELEMENT*>
1301  (base_mesh_pt->element_pt(e));
1302 
1303  // Flush information associated with the external elements
1304  el_pt->flush_all_external_element_storage();
1305 
1306  // Disable projection
1307  el_pt->disable_projection();
1308  }
1309 
1310  //Now unpin everything to restore the problem to its virgin state
1311  this->unpin_all();
1312 
1313  // Now cleanup the storage
1314  Solid_backup.clear();
1315 
1316  /* unsigned ndof_tmp=this->assign_eqn_numbers(); */
1317  if (!Output_during_projection_suppressed)
1318  {
1319  /* oomph_info << "Number of unknowns after project: " */
1320  /* << ndof_tmp << std::endl; */
1321  //std::ostringstream warn_message;
1322  //warn_message
1323  // << "WARNING: Ensure to assign equations numbers in the new mesh,\n"
1324  // << "this is done by calling the assign_eqn_numbers() method from\n"
1325  // << "the original Problem object that has an instance of the mesh.\n"
1326  // << "This may be performed in actions_after_adapt() if the projection\n"
1327  // << "was performed as part of the mesh adaptation process\n\n";
1328  //OomphLibWarning(warn_message.str(),
1329  // OOMPH_CURRENT_FUNCTION,
1330  // OOMPH_EXCEPTION_LOCATION);
1331  }
1332  else
1333  {
1334  // Reset verbosity in Newton solver
1335  Shut_up_in_newton_solve=shut_up_in_newton_solve_backup;
1336 
1337  /// \short Disable documentation of solve times
1338  if (backed_up_doc_time_enabled)
1339  {
1340  linear_solver_pt()->enable_doc_time();
1341  }
1342  }
1343 
1344  } //End of function Projection
1345 
1346  private:
1347 
1348  ///Default constructor (made this private so only the friend class
1349  ///can call the constructor)
1351  {
1352  //This is a linear problem so avoid checking the residual
1353  //after the solve
1354  Problem_is_nonlinear=false; // DO NOT CHANGE THIS -- EVER -- IN
1355  // SOLID MECHANICS PROBLEMS
1356 
1357  // By default suppress output during projection
1358  Output_during_projection_suppressed=true;
1359 
1360  // By default we use an iterative solver for projection
1361  Use_iterative_solver_for_projection=true;
1362 
1363  // Initialise the pointer to the solver and the preconditioner
1364  Iterative_solver_projection_pt = 0;
1365  Preconditioner_projection_pt = 0;
1366  }
1367 
1368  // Destructor
1370  {
1371  if (Iterative_solver_projection_pt!=0)
1372  {
1373  // Clean up memory
1374  delete Iterative_solver_projection_pt;
1375  Iterative_solver_projection_pt = 0;
1376  }
1377 
1378  if (Preconditioner_projection_pt!=0)
1379  {
1380  delete Preconditioner_projection_pt;
1381  Preconditioner_projection_pt = 0;
1382  }
1383 
1384  }
1385 
1386 
1387  /// \short Helper function to store positions (the only things that
1388  /// have been set before doing projection
1390  {
1391  // No need to do anything if there are no elements (in fact, we
1392  // probably never get here...)
1393  if (Problem::mesh_pt()->nelement()==0) return;
1394 
1395  // Deal with positional dofs if (pseudo-)solid element
1396  //If we can cast the first element to a SolidFiniteElement then
1397  //assume that we have a solid mesh
1398  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1399  (Problem::mesh_pt()->element_pt(0));
1400  if(solid_el_pt!=0)
1401  {
1402  const unsigned n_node = this->mesh_pt()->nnode();
1403  Solid_backup.resize(n_node);
1404  //Read dimension and number of position values from the first node
1405  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1406  const unsigned n_position_type =
1407  this->mesh_pt()->node_pt(0)->nposition_type();
1408  //Loop over the nodes
1409  for (unsigned n=0;n<n_node;n++)
1410  {
1411  //Cache a pointer to a solid node
1412  SolidNode* const solid_nod_pt =
1413  dynamic_cast<SolidNode*>(this->mesh_pt()->node_pt(n));
1414  //Now resize the appropriate storage
1415  Solid_backup[n].resize(n_position_type,n_dim);
1416 
1417  for (unsigned i=0;i<n_dim;i++)
1418  {
1419  for(unsigned k=0;k<n_position_type;k++)
1420  {
1421  Solid_backup[n](k,i)= solid_nod_pt->x_gen(k,i);
1422  }
1423  }
1424  }
1425  }
1426  }
1427 
1428  /// \short Helper function to restore positions (the only things that
1429  /// have been set before doing projection
1431  {
1432  // No need to do anything if there are no elements (in fact, we
1433  // probably never get here...)
1434  if (Problem::mesh_pt()->nelement()==0) return;
1435 
1436  // Deal with positional dofs if (pseudo-)solid element
1437  //If we can cast the first element to a SolidFiniteElement then
1438  //assume that we have a solid mesh
1439  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1440  (Problem::mesh_pt()->element_pt(0));
1441  if(solid_el_pt!=0)
1442  {
1443  const unsigned n_node = this->mesh_pt()->nnode();
1444  //Read dimension and number of position values from the first node
1445  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1446  const unsigned n_position_type =
1447  this->mesh_pt()->node_pt(0)->nposition_type();
1448  //Loop over the nodes
1449  for (unsigned n=0;n<n_node;n++)
1450  {
1451  //Cache a pointer to a solid node
1452  SolidNode* const solid_nod_pt =
1453  dynamic_cast<SolidNode*>(this->mesh_pt()->node_pt(n));
1454 
1455  for (unsigned i=0;i<n_dim;i++)
1456  {
1457  for(unsigned k=0;k<n_position_type;k++)
1458  {
1459  solid_nod_pt->x_gen(k,i) = Solid_backup[n](k,i);
1460  }
1461  }
1462  }
1463  }
1464  }
1465 
1466  ///\short Pin all the field values and position unknowns (bit inefficient)
1467  void pin_all()
1468  {
1469  // No need to do anything if there are no elements (in fact, we
1470  // probably never get here...)
1471  if (Problem::mesh_pt()->nelement()==0) return;
1472 
1473  //Loop over all the elements
1474  const unsigned n_element = Problem::mesh_pt()->nelement();
1475  for(unsigned e=0;e<n_element;++e)
1476  {
1477  FiniteElement* el_pt=Problem::mesh_pt()->finite_element_pt(e);
1478  unsigned nint=el_pt->ninternal_data();
1479  for (unsigned j=0;j<nint;j++)
1480  {
1481  Data* int_data_pt = el_pt->internal_data_pt(j);
1482  unsigned nval=int_data_pt->nvalue();
1483  for (unsigned i=0;i<nval;i++)
1484  {
1485  int_data_pt->pin(i);
1486  }
1487  }
1488 
1489  unsigned nnod=el_pt->nnode();
1490  for (unsigned j=0;j<nnod;j++)
1491  {
1492  Node* nod_pt = el_pt->node_pt(j);
1493  unsigned nval=nod_pt->nvalue();
1494  for (unsigned i=0;i<nval;i++)
1495  {
1496  nod_pt->pin(i);
1497  }
1498  }
1499  }
1500 
1501  /// Do we have a solid mesh?
1502  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1503  (Problem::mesh_pt()->element_pt(0));
1504  if (solid_el_pt!=0)
1505  {
1506  //Find number of nodes
1507  const unsigned n_node=this->mesh_pt()->nnode();
1508  //If no nodes then return
1509  if(n_node==0) {return;}
1510 
1511  //Read dimension and number of position values from the first node
1512  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1513  const unsigned n_position_type =
1514  this->mesh_pt()->node_pt(0)->nposition_type();
1515 
1516  //Loop over the nodes
1517  for (unsigned n=0;n<n_node;n++)
1518  {
1519  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(
1520  this->mesh_pt()->node_pt(n));
1521  for (unsigned i=0;i<n_dim;i++)
1522  {
1523  for(unsigned k=0;k<n_position_type;k++)
1524  {
1525  solid_nod_pt->pin_position(k,i);
1526  }
1527  }
1528  }
1529  }
1530  }
1531 
1532 
1533  ///\short Unpin all the field values and position unknowns (bit inefficient)
1534  void unpin_all()
1535  {
1536  // No need to do anything if there are no elements (in fact, we
1537  // probably never get here...)
1538  if (Problem::mesh_pt()->nelement()==0) return;
1539 
1540  //Loop over all the elements
1541  const unsigned n_element = Problem::mesh_pt()->nelement();
1542  for(unsigned e=0;e<n_element;++e)
1543  {
1544  //Cast the first element
1545  PROJECTABLE_ELEMENT * new_el_pt =
1546  dynamic_cast<PROJECTABLE_ELEMENT*>
1547  (Problem::mesh_pt()->element_pt(e));
1548  //Find the number of fields
1549  unsigned n_fields = new_el_pt->nfields_for_projection();
1550 
1551  //Now loop over all fields
1552  for(unsigned f=0;f<n_fields;f++)
1553  {
1554  //Extract the data and value for the field
1555  Vector<std::pair<Data*,unsigned> >
1556  data=new_el_pt->data_values_of_field(f);
1557  //Now loop over all the data and unpin the appropriate values
1558  unsigned d_size=data.size();
1559  for(unsigned d=0;d<d_size;d++)
1560  {
1561  data[d].first->unpin(data[d].second);
1562  }
1563  }
1564  }
1565 
1566  /// Do we have a solid mesh?
1567  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1568  (Problem::mesh_pt()->element_pt(0));
1569  if (solid_el_pt!=0)
1570  {
1571  //Find number of nodes
1572  const unsigned n_node=this->mesh_pt()->nnode();
1573  //If no nodes then return
1574  if(n_node==0) {return;}
1575 
1576  //Read dimension and number of position values from the first node
1577  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1578  const unsigned n_position_type =
1579  this->mesh_pt()->node_pt(0)->nposition_type();
1580 
1581  //Loop over the nodes
1582  for (unsigned n=0;n<n_node;n++)
1583  {
1584  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(
1585  this->mesh_pt()->node_pt(n));
1586  for (unsigned i=0;i<n_dim;i++)
1587  {
1588  for(unsigned k=0;k<n_position_type;k++)
1589  {
1590  solid_nod_pt->unpin_position(k,i);
1591  }
1592  }
1593  }
1594  }
1595  }
1596 
1597  /// \short Helper function to unpin the values of coordinate coord
1598  void unpin_dofs_of_coordinate(const unsigned &coord)
1599  {
1600  //Loop over the nodes
1601  const unsigned n_node = Problem::mesh_pt()->nnode();
1602  //If there are no nodes return immediately
1603  if(n_node==0) {return;}
1604 
1605  //Find the number of position values (should be the same for all nodes)
1606  const unsigned n_position_type =
1607  Problem::mesh_pt()->node_pt(0)->nposition_type();
1608 
1609  for(unsigned n=0;n<n_node;++n)
1610  {
1611  //Cache access to the Node as a solid node
1612  SolidNode* solid_nod_pt =
1613  static_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
1614  //Unpin all position types associated with the given coordinate
1615  for(unsigned k=0;k<n_position_type;++k)
1616  {
1617  solid_nod_pt->unpin_position(k,coord);
1618  }
1619  }
1620  }
1621 
1622  /// \short Helper function to unpin the values of coordinate coord
1623  void pin_dofs_of_coordinate(const unsigned &coord)
1624  {
1625  //Loop over the nodes
1626  const unsigned n_node = Problem::mesh_pt()->nnode();
1627  //If there are no nodes return immediately
1628  if(n_node==0) {return;}
1629 
1630  //Find the number of position values (should be the same for all nodes)
1631  const unsigned n_position_type =
1632  Problem::mesh_pt()->node_pt(0)->nposition_type();
1633 
1634  for(unsigned n=0;n<n_node;++n)
1635  {
1636  //Cache access to the Node as a solid node
1637  SolidNode* solid_nod_pt =
1638  static_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
1639  //Unpin all position types associated with the given coordinate
1640  for(unsigned k=0;k<n_position_type;++k)
1641  {
1642  solid_nod_pt->pin_position(k,coord);
1643  }
1644  }
1645  }
1646 
1647 
1648  ///Helper function to unpin dofs of fld-th field
1649  void unpin_dofs_of_field(const unsigned &fld)
1650  {
1651  unsigned n_element = Problem::mesh_pt()->nelement();
1652  for(unsigned e=0;e<n_element;e++)
1653  {
1654  PROJECTABLE_ELEMENT * new_el_pt =
1655  dynamic_cast<PROJECTABLE_ELEMENT*>
1656  (Problem::mesh_pt()->element_pt(e));
1657 
1658  Vector<std::pair<Data*,unsigned> >
1659  data=new_el_pt->data_values_of_field(fld);
1660 
1661  unsigned d_size=data.size();
1662  for(unsigned d=0;d<d_size;d++)
1663  {
1664  data[d].first->unpin(data[d].second);
1665  }
1666  }
1667  }
1668 
1669  ///Helper function to unpin dofs of fld-th field
1670  void pin_dofs_of_field(const unsigned &fld)
1671  {
1672  unsigned n_element = Problem::mesh_pt()->nelement();
1673  for(unsigned e=0;e<n_element;e++)
1674  {
1675  PROJECTABLE_ELEMENT * new_el_pt =
1676  dynamic_cast<PROJECTABLE_ELEMENT*>
1677  (Problem::mesh_pt()->element_pt(e));
1678 
1679  Vector<std::pair<Data*,unsigned> >
1680  data=new_el_pt->data_values_of_field(fld);
1681 
1682  unsigned d_size=data.size();
1683  for(unsigned d=0;d<d_size;d++)
1684  {
1685  data[d].first->pin(data[d].second);
1686  }
1687  }
1688  }
1689 
1690  /// Helper function to set time level for projection
1691  void set_time_level_for_projection(const unsigned &time_level)
1692  {
1693  unsigned n_element = Problem::mesh_pt()->nelement();
1694  for(unsigned e=0;e<n_element;e++)
1695  {
1696  PROJECTABLE_ELEMENT * el_pt =
1697  dynamic_cast<PROJECTABLE_ELEMENT*>
1698  (Problem::mesh_pt()->element_pt(e));
1699 
1700  //Set what time we are dealing with
1701  el_pt->time_level_for_projection()=time_level;
1702  }
1703  }
1704 
1705  ///Set the coordinate for projection
1706  void set_coordinate_for_projection(const unsigned &i)
1707  {
1708  unsigned n_element = Problem::mesh_pt()->nelement();
1709  for(unsigned e=0;e<n_element;e++)
1710  {
1711  PROJECTABLE_ELEMENT * new_el_pt =
1712  dynamic_cast<PROJECTABLE_ELEMENT*>
1713  (Problem::mesh_pt()->element_pt(e));
1714 
1715  //Set that we are solving a projected coordinate problem
1716  new_el_pt->set_project_coordinates();
1717 
1718  //Set the projected coordinate
1719  new_el_pt->projected_coordinate()=i;
1720  }
1721  }
1722 
1723  ///Set the Lagrangian coordinate for projection
1725  {
1726  unsigned n_element = Problem::mesh_pt()->nelement();
1727  for(unsigned e=0;e<n_element;e++)
1728  {
1729  PROJECTABLE_ELEMENT * new_el_pt =
1730  dynamic_cast<PROJECTABLE_ELEMENT*>
1731  (Problem::mesh_pt()->element_pt(e));
1732 
1733  //Set that we are solving a projected Lagrangian coordinate problem
1734  new_el_pt->set_project_lagrangian();
1735 
1736  //Set the projected coordinate
1737  new_el_pt->projected_lagrangian_coordinate()=i;
1738  }
1739  }
1740 
1741  ///Set current field for projection
1742  void set_current_field_for_projection(const unsigned &fld)
1743  {
1744  unsigned n_element = Problem::mesh_pt()->nelement();
1745  for(unsigned e=0;e<n_element;e++)
1746  {
1747  PROJECTABLE_ELEMENT * new_el_pt =
1748  dynamic_cast<PROJECTABLE_ELEMENT*>
1749  (Problem::mesh_pt()->element_pt(e));
1750 
1751  //Set current field
1752  new_el_pt->projected_field()=fld;
1753  }
1754  }
1755 
1756  private:
1757 
1758  /// Backup for pin status of solid node's position Data
1759  Vector<DenseMatrix<double> > Solid_backup;
1760 
1761  /// Flag to suppress output during projection
1763 
1764  // Use an iterative solver for solving the system of equations
1766 
1767  // The iterative solver to solve the projection problem
1768  IterativeLinearSolver *Iterative_solver_projection_pt;
1769 
1770  // The preconditioner for the solver
1772 
1773 };
1774 
1775 
1776 } //End of namespace
1777 
1778 #endif
unsigned Backup_ninteraction
Store number of "external" interactions that were assigned to the element before doing the projection...
Definition: projection.h:95
Projection_Type Projection_type
Variable to indicate if we&#39;re projecting the history values of the nodal coordinates (Coordinate) the...
Definition: projection.h:85
void set_coordinate_for_projection(const unsigned &i)
Set the coordinate for projection.
Definition: projection.h:1706
void disable_projection()
Helper function to restore the element to the state it was in before we entered the projection mode a...
Definition: projection.h:583
void unpin_dofs_of_coordinate(const unsigned &coord)
Helper function to unpin the values of coordinate coord.
Definition: projection.h:1598
unsigned & time_level_for_projection()
Which history value are we projecting?
Definition: projection.h:627
IterativeLinearSolver * Iterative_solver_projection_pt
Definition: projection.h:1768
void pin_all()
Pin all the field values and position unknowns (bit inefficient)
Definition: projection.h:1467
ProjectableElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
Definition: projection.h:238
virtual double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)=0
Return Jacobian of mapping and the shape functions associated with field fld. The number of shape fun...
void set_project_lagrangian()
Project (current and only values of) lagrangian coordinates.
Definition: projection.h:619
bool Do_projection
Bool to know if we do projection or not. If false (the default) we solve the element&#39;s "real" equatio...
Definition: projection.h:90
void pin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
Definition: projection.h:1670
Template-free Base class for projectable elements.
Definition: projection.h:59
void enable_suppress_output_during_projection()
Suppress all output during projection phases.
Definition: projection.h:701
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Overloaded version of fill_in_contribution_to_residuals.
Definition: projection.h:187
unsigned & projected_field()
Field that is currently being projected.
Definition: projection.h:623
void set_project_coordinates()
Project (history values of) coordintes.
Definition: projection.h:613
void project(Mesh *base_mesh_pt, const bool &dont_project_positions=false)
Project from base into the problem&#39;s own mesh.
Definition: projection.h:726
unsigned Time_level_for_projection
Time level we are projecting (0=current values; >0: history values)
Definition: projection.h:72
virtual ~ProjectableElementBase()
Virtual destructor.
Definition: projection.h:120
void set_lagrangian_coordinate_for_projection(const unsigned &i)
Set the Lagrangian coordinate for projection.
Definition: projection.h:1724
void store_positions()
Helper function to store positions (the only things that have been set before doing projection...
Definition: projection.h:1389
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we&#39;re projecting...
Definition: projection.h:76
bool Output_during_projection_suppressed
Flag to suppress output during projection.
Definition: projection.h:1762
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Overloaded version of fill_in_contribution_to_jacobian.
Definition: projection.h:218
void unpin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
Definition: projection.h:1649
bool use_iterative_solver_for_projection()
Definition: projection.h:714
virtual void residual_for_projection(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Residual for the projection step. Flag indicates if we want the Jacobian (1) or not (0)...
Definition: projection.h:243
virtual unsigned nhistory_values_for_projection(const unsigned &fld)=0
Number of history values to be stored for fld-th field (includes current value!)
virtual unsigned nvalue_of_field(const unsigned &fld)=0
Return number of values (pinned or not) associated with field fld within the element. This must correspond to the number of shape functions returned in jacobian_and_shape_of_field(...).
virtual unsigned nhistory_values_for_coordinate_projection()=0
Number of history values to be stored when projecting the history values of the nodal coordinates (in...
void enable_projection()
Backup the element&#39;s state and switch it to projection mode.
Definition: projection.h:548
bool Backup_external_interaction_data
Remember if the element includes external data when used in non-projection mode (this is temporarily ...
Definition: projection.h:106
void enable_use_iterative_solver_for_projection()
Enables the use of an iterative solver for projection.
Definition: projection.h:718
unsigned & projected_coordinate()
When projecting the history values of the nodal coordinates, this is the coordinate we&#39;re projecting...
Definition: projection.h:632
void pin_dofs_of_coordinate(const unsigned &coord)
Helper function to unpin the values of coordinate coord.
Definition: projection.h:1623
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Use Eulerian coordinates for matching in locate_zeta when doing projection.
Definition: projection.h:531
unsigned & projected_lagrangian_coordinate()
When projecting the Lagrangian coordinates this is the coordinate that is being projected.
Definition: projection.h:637
void unpin_all()
Unpin all the field values and position unknowns (bit inefficient)
Definition: projection.h:1534
virtual unsigned nfields_for_projection()=0
Number of fields of the problem, so e.g. for 2D Navier Stokes this would be 3 (for the two velocities...
void set_time_level_for_projection(const unsigned &time_level)
Helper function to set time level for projection.
Definition: projection.h:1691
void disable_use_iterative_solver_for_projection()
Disbales the use of an iterative solver for projection.
Definition: projection.h:722
void set_current_field_for_projection(const unsigned &fld)
Set current field for projection.
Definition: projection.h:1742
Vector< DenseMatrix< double > > Solid_backup
Backup for pin status of solid node&#39;s position Data.
Definition: projection.h:1759
virtual double get_field(const unsigned &time, const unsigned &fld, const Vector< double > &s)=0
Return the fld-th field at local coordinates s at time-level time (time=0: current value; time>0: his...
void disable_suppress_output_during_projection()
Undo suppression of all output during projection phases.
Definition: projection.h:707
void restore_positions()
Helper function to restore positions (the only things that have been set before doing projection...
Definition: projection.h:1430
virtual Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)=0
Pure virtual function in which the element writer must specify the values associated with field fld...
unsigned Projected_field
Field that is currently being projected.
Definition: projection.h:69
virtual int local_equation(const unsigned &fld, const unsigned &ivalue)=0
Return local equation numbers associated with value ivalue of field fld within the element...
bool Use_iterative_solver_for_projection
Definition: projection.h:1765
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Definition: projection.h:210
Projection_Type
Enumerated collection to specify which projection problem is to be solved.
Definition: projection.h:66
bool Backup_external_geometric_data
Remember if the element includes external geometric data when used in non-projection mode (this is te...
Definition: projection.h:100
Preconditioner * Preconditioner_projection_pt
Definition: projection.h:1771
void set_project_values()
Project (history values of) values.
Definition: projection.h:616
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.
Definition: projection.h:80