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