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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
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"
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
48 
49 // Use a preconditioner for the iterative solver
50 #include "preconditioner.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.
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.
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
188  {
189  //Do projection
190  if (Do_projection)
191  {
193  (residuals,GeneralisedElement::Dummy_matrix, 0);
194  }
195  //solve problem normally
196  else
197  {
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  {
214  ELEMENT::describe_local_dofs(out,current_string);
215  }
216 
217  ///Overloaded version of fill_in_contribution_to_jacobian
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=
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
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
552 
553  //Backup flag for inclusion of geometric data
555  {
557  }
558  else
559  {
561  }
562 
563  //Backup flag for inclusion of interaction data
565  {
567  }
568  else
569  {
571  }
572 
573  //Actions to enable projection
574  Do_projection=true;
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
587 
588  //Restore geometric data
590  {
592  }
593  else
594  {
596  }
597 
598  //Restore interaction data
600  {
602  }
603  else
604  {
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>
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  // The classes are friend only when the templated element of the
690  // friend class matches the templated element of the
691  // ProjectionProblem class
692  // friend class RefineableTriangleMesh<class PROJECTABLE_ELEMENT>;
693  // friend class RefineableTetgenMesh<class PROJECTABLE_ELEMENT>;
694  // friend class BackupMeshForProjection<class PROJECTABLE_ELEMENT>;
695 
696  public:
697 
698  /// Suppress all output during projection phases
700  {
702  }
703 
704  /// Undo suppression of all output during projection phases
706  {
708  }
709 
710  /// Return the value of the flag about using an iterative solver for
711  /// projection
714 
715  /// Enables the use of an iterative solver for projection
718 
719  /// Disbales the use of an iterative solver for projection
722 
723  ///\short Project from base into the problem's own mesh.
724  void project(Mesh* base_mesh_pt, const bool& dont_project_positions=false)
725  {
726  // Use an iterative solver?
728  {
729  // If oomph-lib has Trilinos installed then use the CG version
730  // from Trilinos, otherwise use oomph-lib's own CG implementation
731 #ifdef OOMPH_HAS_TRILINOS
732  // Check whether the problem is distributed?
734  {
735  // Create a Trilinos Solver
737  // Select CG as the linear solver
739  solver_type() = TrilinosAztecOOSolver::CG;
740  }
741  else
742  {
743  // Use CG to solve the projection problem
745  }
746 
747  // Create the preconditioner
749  // Set the preconditioner for the solver
752 
753  // Set CG as the linear solver
755 
756 #else
757  // Check whether the problem is distributed?
759  {
760  // If we did not installed Trilinos and the problem is not
761  // distributed then we can use a (serial) preconditioned
762  // iterative solver, otherwise, if we did not installed Trilinos
763  // but the problem is distributed then we cannot use a
764  // preconditioned iterative solver. Matrix multiplication in a
765  // distributed environment is only performed by Trilinos. We
766  // then use a direct solver for the projection problem.
767 
768  // Use CG to solve the projection problem
770 
771  // Create the preconditioner
773  // Set the preconditioner for the solver
776 
777  // Set CG as the linear solver
779  }
780  else
781  {
782  // Use a direct solver. Do nothing
783  }
784 
785 #endif
786 
787  } // if (Use_iterative_solver_for_projection)
788 
789  // Backup verbosity in Newton solve status
790  bool shut_up_in_newton_solve_backup=Shut_up_in_newton_solve;
791 
792  // Disable documentation of solve times
793  bool backed_up_doc_time_enabled=linear_solver_pt()->is_doc_time_enabled();
795  {
797  }
798 
799  //Display stats
800  unsigned n_element = Problem::mesh_pt()->nelement();
801  unsigned n_element1=base_mesh_pt->nelement();
802  unsigned n_node = Problem::mesh_pt()->nnode();
803  unsigned n_node1=base_mesh_pt->nnode();
805  {
806  oomph_info <<"\n=============================\n";
807  oomph_info << "Base mesh has " << n_element1 << " elements\n";
808  oomph_info << "Target mesh has " << n_element << " elements\n";
809  oomph_info << "Base mesh has " << n_node1 << " nodes\n";
810  oomph_info << "Target mesh has " << n_node << " nodes\n";
811  oomph_info <<"=============================\n\n";
812 
813  }
814  else
815  {
816  // Make Newton solver shut up too
818  }
819 
820 
821  if (n_element==0)
822  {
823  oomph_info
824  << "Very odd -- no elements in target mesh; "
825  << " not doing anything in ProjectionProblem::project()\n";
826  return;
827  }
828 
829 #ifdef PARANOID
830  unsigned nnod=Problem::mesh_pt()->nnode();
831  if (nnod==0)
832  {
833  std::ostringstream error_stream;
834  error_stream
835  << "Mesh has no nodes! Please populate the Node_pt vector\n"
836  << "otherwise projection won't work!\n";
837  throw OomphLibError(error_stream.str(),
838  OOMPH_CURRENT_FUNCTION,
839  OOMPH_EXCEPTION_LOCATION);
840  }
841 #endif
842 
843  // How many fields do we have to project?
844  unsigned n_fields = dynamic_cast<PROJECTABLE_ELEMENT*>
845  (Problem::mesh_pt()->element_pt(0))->nfields_for_projection();
846 
847  // Spatial dimension of the problem
848  unsigned n_dim = Problem::mesh_pt()->node_pt(0)->ndim();
849 
850  // Default number of history values
851  unsigned n_history_values=0;
852 
853  // Set things up for coordinate projection
854  for(unsigned e=0;e<n_element;e++)
855  {
856  PROJECTABLE_ELEMENT * el_pt =
857  dynamic_cast<PROJECTABLE_ELEMENT*>
859 
860  // Switch to projection
861  el_pt->enable_projection();
862  }
863 
864  // Switch elements in base mesh to projection mode (required
865  // to switch to use of Eulerian coordinates when identifying
866  // corresponding points in the two meshes)
867  for(unsigned e=0;e<n_element1;e++)
868  {
869  PROJECTABLE_ELEMENT * el_pt =
870  dynamic_cast<PROJECTABLE_ELEMENT*>
871  (base_mesh_pt->element_pt(e));
872 
873  // Switch to projection
874  el_pt->enable_projection();
875  }
876 
877 
878  // Set up multi domain interactions so we can locate the
879  // values in the base mesh.
880  // Note that it's important to switch elements to projection
881  // mode first to ensure that matching is done based on Eulerian
882  // rather than Lagrangian coordinates if pseudo-solid elements
883  // are used.
884  double t_start=TimingHelpers::timer();
886  <PROJECTABLE_ELEMENT>(this,Problem::mesh_pt(),base_mesh_pt);
888  {
889  oomph_info <<"CPU for setup of multi-domain interaction for projection: "
890  << TimingHelpers::timer()-t_start << std::endl;
891  }
892  t_start=TimingHelpers::timer();
893 
894 
895  //Let us first pin every degree of freedom
896  //We shall unpin selected dofs for each different projection problem
897  this->pin_all();
898 
899  if (!dont_project_positions)
900  {
901 
902  //------------------Project coordinates first------------------------
903  //If we have a solid element then we should also project Lagrangian
904  //coordinates, but we can use the storage that MUST be provided for
905  //the unknown positions for this.
906  //If we can cast the first element of the mesh to a solid element,
907  //then assume that we have solid elements
908  if(dynamic_cast<SolidFiniteElement*>(
909  Problem::mesh_pt()->element_pt(0)))
910  {
911  // Store current positions
912  this->store_positions();
913 
914  //There are no history values for the Lagrangian coordinates
915  //Set coordinate 0 for projection
917  this->unpin_dofs_of_coordinate(0);
918 
919  //Loop over the Lagrangian coordinate
920  const unsigned n_lagrangian =
921  dynamic_cast<SolidNode*>(Problem::mesh_pt()->node_pt(0))->nlagrangian();
922 
923  //Now loop over the lagrangian coordinates
924  for(unsigned i=0;i<n_lagrangian;++i)
925  {
926 
928  {
929  oomph_info <<"\n\n=============================================\n";
930  oomph_info << "Projecting values for Lagrangian coordinate " << i
931  << std::endl;
932  oomph_info <<"=============================================\n\n";
933  }
934 
935  //Set the coordinate for projection
937 
938  //Assign equation number
939  unsigned ndof_tmp=assign_eqn_numbers();
941  {
942  oomph_info <<
943  "Number of equations for projection of Lagrangian coordinate "
944  << " : "<< ndof_tmp <<std::endl << std::endl;
945  }
946 
947 
949  {
950  std::ostringstream error_stream;
951  error_stream
952  << "Solid mechanics problems will break if Problem_is_nonlinear is\n"
953  << "set to true in the projection problem because we're using the\n "
954  << "actual nodal positions to store the values of the Lagrangian\n"
955  << "coords. There shouldn't be any need for \n"
956  << "Problem_is_nonlinear=true anyway, apart from debugging in \n"
957  << "which case you now know why this case will break!\n";
958  OomphLibWarning(error_stream.str(),
959  OOMPH_CURRENT_FUNCTION,
960  OOMPH_EXCEPTION_LOCATION);
961  }
962 
963 
964  //Projection and interpolation
966 
967  //Move values back into Lagrangian coordinate for all nodes
968  unsigned n_node = Problem::mesh_pt()->nnode();
969  for(unsigned n=0;n<n_node;++n)
970  {
971  //Cast it to a solid node
972  SolidNode* solid_node_pt =
973  dynamic_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
974  //Now find number of types
975  const unsigned n_position_type = solid_node_pt->nposition_type();
976  //Find number of lagrangian types
977  const unsigned n_lagrangian_type = solid_node_pt->nlagrangian_type();
978 
979  //If these are not the same, throw an error
980  if(n_position_type != n_lagrangian_type)
981  {
982  std::ostringstream error_stream;
983  error_stream
984  << "The number of generalised position dofs "
985  << n_position_type
986  << "\n not the same as the number of generalised lagrangian dofs "
987  << n_lagrangian_type << ".\n"
988  << "Projection cannot be done at the moment, sorry.\n";
989 
990  throw OomphLibError(error_stream.str(),
991  OOMPH_CURRENT_FUNCTION,
992  OOMPH_EXCEPTION_LOCATION);
993  }
994 
995  //Now transfer the information across
996  //from the first coordinate which was used during the projection
997  for(unsigned k=0;k<n_position_type;++k)
998  {
999  solid_node_pt->xi_gen(k,i) = solid_node_pt->x_gen(k,0);
1000  //Reset real values so that the Jacobians are correctly
1001  //computed next time around
1002  solid_node_pt->x_gen(k,0) = Solid_backup[n](k,0);
1003  }
1004  }
1005  } //End of loop over lagrangian coordinates
1006 
1007  //Now repin the dofs
1008  this->pin_dofs_of_coordinate(0);
1009 
1010  //Now project the position histories
1011 
1012  //Check number of history values for coordinates
1013  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1014  (Problem::mesh_pt()->element_pt(0))->
1015  nhistory_values_for_coordinate_projection();
1016 
1017  //Projection the coordinates only if there are history values
1018  if(n_history_values > 1)
1019  {
1020  for (unsigned i=0;i<n_dim;i++)
1021  {
1023  {
1024  oomph_info <<"\n\n=============================================\n";
1025  oomph_info << "Projecting history values for coordinate " << i
1026  << std::endl;
1027  oomph_info <<"=============================================\n\n";
1028  }
1029 
1030  //Set the coordinate for projection
1032  this->unpin_dofs_of_coordinate(i);
1033 
1034  // Loop over number of history values, beginning with the latest one.
1035  // Don't deal with current time.
1036  for (unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1037  {
1038  unsigned time_level=h_tim-1;
1039 
1040  //Set time_level we are dealing with
1041  this->set_time_level_for_projection(time_level);
1042 
1043  //Assign equation number
1044  unsigned ndof_tmp=assign_eqn_numbers();
1046  {
1047  oomph_info << "Number of equations for projection of coordinate "
1048  << i << " at time level " << time_level
1049  << " : "<< ndof_tmp <<std::endl << std::endl;
1050  }
1051 
1052  //Projection and interpolation
1054 
1055  //Move values back into history value of coordinate
1056  unsigned n_node = Problem::mesh_pt()->nnode();
1057  for(unsigned n=0;n<n_node;++n)
1058  {
1059  //Cache the pointer to the node
1060  Node* nod_pt = Problem::mesh_pt()->node_pt(n);
1061  //Find the number of generalised dofs
1062  const unsigned n_position_type = nod_pt->nposition_type();
1063  //Now copy all back
1064  for(unsigned k=0;k<n_position_type;++k)
1065  {
1066  nod_pt->x_gen(time_level,k,i) = nod_pt->x_gen(0,k,i);
1067  //Reset real values so that the Jacobians are computed
1068  //correctly next time around
1069  nod_pt->x_gen(0,k,i) = Solid_backup[n](k,i);
1070  }
1071  }
1072  }
1073  //Repin the positions
1074  this->pin_dofs_of_coordinate(i);
1075  }
1076  } //End of history value projection
1077  } //End of SolidElement case
1078 
1079  //Now for non solid elements, we are going to hijack the
1080  //first value as storage to be used for the projection of the history
1081  //values
1082  else
1083  {
1084  // Prepare for projection in value 0
1086  this->unpin_dofs_of_field(0);
1087 
1088 #ifdef PARANOID
1089 
1090  // The machinery used below assumes that field 0 can actually
1091  // hold the coordinates i.e. that the field is interpolated
1092  // isoparametrically! The method will fail if there are no values
1093  // stored at the nodes. Currently there are no examples of that -- the
1094  // problem would only arise for elements that have all their fields
1095  // represented by internal data. Will fix this if/when
1096  // such a case arises...
1097  unsigned n_element = Problem::mesh_pt()->nelement();
1098  for(unsigned e=0;e<n_element;e++)
1099  {
1101  unsigned nnod=el_pt->nnode();
1102  for (unsigned j=0;j<nnod;j++)
1103  {
1104  Node* nod_pt=el_pt->node_pt(j);
1105  if (nod_pt->nvalue()==0)
1106  {
1107  std::ostringstream error_stream;
1108  error_stream << "Node at ";
1109  unsigned n=nod_pt->ndim();
1110  for (unsigned i=0;i<n;i++)
1111  {
1112  error_stream << nod_pt->x(i) << " ";
1113  }
1114  error_stream
1115  << "\nhas no values. Projection (of coordinates) doesn't work\n"
1116  << "for such cases (at the moment), sorry! Send us the code\n"
1117  << "where the problem arises and we'll try to implement this\n"
1118  << "(up to now unnecessary) capability...\n";
1119  throw OomphLibError(error_stream.str(),
1120  OOMPH_CURRENT_FUNCTION,
1121  OOMPH_EXCEPTION_LOCATION);
1122  }
1123  }
1124  }
1125 
1126 #endif
1127 
1128  //Check number of history values for coordinates
1129  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1130  (Problem::mesh_pt()->element_pt(0))->
1131  nhistory_values_for_coordinate_projection();
1132 
1133  //Projection the coordinates only if there are history values
1134  if(n_history_values > 1)
1135  {
1136  for (unsigned i=0;i<n_dim;i++)
1137  {
1139  {
1140  oomph_info <<"\n\n=============================================\n";
1141  oomph_info << "Projecting history values for coordinate " << i
1142  << std::endl;
1143  oomph_info <<"=============================================\n\n";
1144  }
1145 
1146  //Set the coordinate for projection
1148 
1149  // Loop over number of history values, beginning with the latest one.
1150  // Don't deal with current time.
1151  for (unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1152  {
1153  unsigned time_level=h_tim-1;
1154 
1155  //Set time_level we are dealing with
1156  this->set_time_level_for_projection(time_level);
1157 
1158  //Assign equation number
1159  unsigned ndof_tmp=assign_eqn_numbers();
1161  {
1162  oomph_info << "Number of equations for projection of coordinate "
1163  << i << " at time level " << time_level
1164  << " : "<< ndof_tmp <<std::endl << std::endl;
1165  }
1166 
1167  //Projection and interpolation
1169 
1170  //Move values back into history value of coordinate
1171  unsigned n_element = Problem::mesh_pt()->nelement();
1172  for(unsigned e=0;e<n_element;e++)
1173  {
1174  PROJECTABLE_ELEMENT * new_el_pt =
1175  dynamic_cast<PROJECTABLE_ELEMENT*>
1177 
1179  data=new_el_pt->data_values_of_field(0);
1180 
1181  unsigned d_size=data.size();
1182  for(unsigned d=0;d<d_size;d++)
1183  {
1184  //Replace as coordinates
1185  double coord=data[d].first->value(0,0);
1186  dynamic_cast<Node*>(data[d].first)->x(time_level,i) = coord;
1187  }
1188  }
1189  }
1190  }
1191  } //End of history value projection
1192 
1193  //Repin the dofs for field 0
1194  this->pin_dofs_of_field(0);
1195 
1196  } //End of non-SolidElement case
1197 
1198 
1199  } // end if for projection of coordinates
1200 
1201  //Disable projection of coordinates
1202  for(unsigned e=0;e<n_element;e++)
1203  {
1204  PROJECTABLE_ELEMENT * el_pt =
1205  dynamic_cast<PROJECTABLE_ELEMENT*>
1207 
1208  el_pt->set_project_values();
1209  }
1210 
1211  //Loop over fields
1212  for (unsigned fld=0; fld<n_fields ;fld++)
1213  {
1214  //Let us first pin every degree of freedom
1215  //We shall unpin selected dofs for each different projection problem
1216  this->pin_all();
1217 
1218  //Do actions for this field
1220  this->unpin_dofs_of_field(fld);
1221 
1222  //Check number of history values
1223  n_history_values = dynamic_cast<PROJECTABLE_ELEMENT*>
1224  (Problem::mesh_pt()->element_pt(0))->nhistory_values_for_projection(fld);
1225 
1226  //Loop over number of history values
1227  //Beginning with the latest one
1228  for (unsigned h_tim=n_history_values;h_tim>0;h_tim--)
1229  {
1230  unsigned time_level=h_tim-1;
1232  {
1233  oomph_info <<"\n=========================================\n";
1234  oomph_info << "Projecting field " << fld << " at time level "
1235  << time_level<<std::endl;
1236  oomph_info << "========================================\n";
1237  }
1238 
1239  //Set time_level we are dealing with
1240  this->set_time_level_for_projection(time_level);
1241 
1242  //Assign equation number
1243  unsigned ndof_tmp=assign_eqn_numbers();
1245  {
1246  oomph_info << "Number of equations for projection of field "
1247  << fld << " at time level " << time_level
1248  << " : "<< ndof_tmp <<std::endl << std::endl;
1249  }
1250 
1251  //Projection and interpolation
1253 
1254  // Move computed values into the required time-level (not needed
1255  // for current values which are done last -- they simply
1256  // stay where they are)
1257  if(time_level!=0)
1258  {
1259  for(unsigned e=0;e<n_element;e++)
1260  {
1261  PROJECTABLE_ELEMENT * new_el_pt =
1262  dynamic_cast<PROJECTABLE_ELEMENT*>
1264 
1266  data=new_el_pt->data_values_of_field(fld);
1267 
1268  unsigned d_size=data.size();
1269  for(unsigned d=0;d<d_size;d++)
1270  {
1271  //Move into time level
1272  double c_value=data[d].first->value(0,data[d].second);
1273  data[d].first->set_value(time_level,data[d].second,c_value);
1274  }
1275  }
1276  }
1277  } //End of loop over time levels
1278 
1279  } //End of loop over fields
1280 
1281 
1282  //Reset parameters of external storage and interactions
1283  for(unsigned e=0;e<n_element;e++)
1284  {
1285  PROJECTABLE_ELEMENT * new_el_pt =
1286  dynamic_cast<PROJECTABLE_ELEMENT*>
1288 
1289  // Flush information associated with the external elements
1290  new_el_pt->flush_all_external_element_storage();
1291 
1292  new_el_pt->disable_projection();
1293  }
1294 
1295  for(unsigned e=0;e<n_element1;e++)
1296  {
1297  PROJECTABLE_ELEMENT * el_pt =
1298  dynamic_cast<PROJECTABLE_ELEMENT*>
1299  (base_mesh_pt->element_pt(e));
1300 
1301  // Flush information associated with the external elements
1302  el_pt->flush_all_external_element_storage();
1303 
1304  // Disable projection
1305  el_pt->disable_projection();
1306  }
1307 
1308  //Now unpin everything to restore the problem to its virgin state
1309  this->unpin_all();
1310 
1311  // Now cleanup the storage
1312  Solid_backup.clear();
1313 
1314  /* unsigned ndof_tmp=this->assign_eqn_numbers(); */
1316  {
1317  /* oomph_info << "Number of unknowns after project: " */
1318  /* << ndof_tmp << std::endl; */
1319  //std::ostringstream warn_message;
1320  //warn_message
1321  // << "WARNING: Ensure to assign equations numbers in the new mesh,\n"
1322  // << "this is done by calling the assign_eqn_numbers() method from\n"
1323  // << "the original Problem object that has an instance of the mesh.\n"
1324  // << "This may be performed in actions_after_adapt() if the projection\n"
1325  // << "was performed as part of the mesh adaptation process\n\n";
1326  //OomphLibWarning(warn_message.str(),
1327  // OOMPH_CURRENT_FUNCTION,
1328  // OOMPH_EXCEPTION_LOCATION);
1329  }
1330  else
1331  {
1332  // Reset verbosity in Newton solver
1333  Shut_up_in_newton_solve=shut_up_in_newton_solve_backup;
1334 
1335  /// \short Disable documentation of solve times
1336  if (backed_up_doc_time_enabled)
1337  {
1339  }
1340  }
1341 
1342  } //End of function Projection
1343 
1344  private:
1345 
1346  ///Default constructor (made this private so only the friend class
1347  ///can call the constructor)
1349  {
1350  //This is a linear problem so avoid checking the residual
1351  //after the solve
1352  Problem_is_nonlinear=false; // DO NOT CHANGE THIS -- EVER -- IN
1353  // SOLID MECHANICS PROBLEMS
1354 
1355  // By default suppress output during projection
1357 
1358  // By default we use an iterative solver for projection
1360 
1361  // Initialise the pointer to the solver and the preconditioner
1364  }
1365 
1366  // Destructor
1368  {
1370  {
1371  // Clean up memory
1374  }
1375 
1377  {
1380  }
1381 
1382  }
1383 
1384 
1385  /// \short Helper function to store positions (the only things that
1386  /// have been set before doing projection
1388  {
1389  // No need to do anything if there are no elements (in fact, we
1390  // probably never get here...)
1391  if (Problem::mesh_pt()->nelement()==0) return;
1392 
1393  // Deal with positional dofs if (pseudo-)solid element
1394  //If we can cast the first element to a SolidFiniteElement then
1395  //assume that we have a solid mesh
1396  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1397  (Problem::mesh_pt()->element_pt(0));
1398  if(solid_el_pt!=0)
1399  {
1400  const unsigned n_node = this->mesh_pt()->nnode();
1401  Solid_backup.resize(n_node);
1402  //Read dimension and number of position values from the first node
1403  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1404  const unsigned n_position_type =
1405  this->mesh_pt()->node_pt(0)->nposition_type();
1406  //Loop over the nodes
1407  for (unsigned n=0;n<n_node;n++)
1408  {
1409  //Cache a pointer to a solid node
1410  SolidNode* const solid_nod_pt =
1411  dynamic_cast<SolidNode*>(this->mesh_pt()->node_pt(n));
1412  //Now resize the appropriate storage
1413  Solid_backup[n].resize(n_position_type,n_dim);
1414 
1415  for (unsigned i=0;i<n_dim;i++)
1416  {
1417  for(unsigned k=0;k<n_position_type;k++)
1418  {
1419  Solid_backup[n](k,i)= solid_nod_pt->x_gen(k,i);
1420  }
1421  }
1422  }
1423  }
1424  }
1425 
1426  /// \short Helper function to restore positions (the only things that
1427  /// have been set before doing projection
1429  {
1430  // No need to do anything if there are no elements (in fact, we
1431  // probably never get here...)
1432  if (Problem::mesh_pt()->nelement()==0) return;
1433 
1434  // Deal with positional dofs if (pseudo-)solid element
1435  //If we can cast the first element to a SolidFiniteElement then
1436  //assume that we have a solid mesh
1437  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1438  (Problem::mesh_pt()->element_pt(0));
1439  if(solid_el_pt!=0)
1440  {
1441  const unsigned n_node = this->mesh_pt()->nnode();
1442  //Read dimension and number of position values from the first node
1443  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1444  const unsigned n_position_type =
1445  this->mesh_pt()->node_pt(0)->nposition_type();
1446  //Loop over the nodes
1447  for (unsigned n=0;n<n_node;n++)
1448  {
1449  //Cache a pointer to a solid node
1450  SolidNode* const solid_nod_pt =
1451  dynamic_cast<SolidNode*>(this->mesh_pt()->node_pt(n));
1452 
1453  for (unsigned i=0;i<n_dim;i++)
1454  {
1455  for(unsigned k=0;k<n_position_type;k++)
1456  {
1457  solid_nod_pt->x_gen(k,i) = Solid_backup[n](k,i);
1458  }
1459  }
1460  }
1461  }
1462  }
1463 
1464  ///\short Pin all the field values and position unknowns (bit inefficient)
1465  void pin_all()
1466  {
1467  // No need to do anything if there are no elements (in fact, we
1468  // probably never get here...)
1469  if (Problem::mesh_pt()->nelement()==0) return;
1470 
1471  //Loop over all the elements
1472  const unsigned n_element = Problem::mesh_pt()->nelement();
1473  for(unsigned e=0;e<n_element;++e)
1474  {
1476  unsigned nint=el_pt->ninternal_data();
1477  for (unsigned j=0;j<nint;j++)
1478  {
1479  Data* int_data_pt = el_pt->internal_data_pt(j);
1480  unsigned nval=int_data_pt->nvalue();
1481  for (unsigned i=0;i<nval;i++)
1482  {
1483  int_data_pt->pin(i);
1484  }
1485  }
1486 
1487  unsigned nnod=el_pt->nnode();
1488  for (unsigned j=0;j<nnod;j++)
1489  {
1490  Node* nod_pt = el_pt->node_pt(j);
1491  unsigned nval=nod_pt->nvalue();
1492  for (unsigned i=0;i<nval;i++)
1493  {
1494  nod_pt->pin(i);
1495  }
1496  }
1497  }
1498 
1499  /// Do we have a solid mesh?
1500  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1501  (Problem::mesh_pt()->element_pt(0));
1502  if (solid_el_pt!=0)
1503  {
1504  //Find number of nodes
1505  const unsigned n_node=this->mesh_pt()->nnode();
1506  //If no nodes then return
1507  if(n_node==0) {return;}
1508 
1509  //Read dimension and number of position values from the first node
1510  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1511  const unsigned n_position_type =
1512  this->mesh_pt()->node_pt(0)->nposition_type();
1513 
1514  //Loop over the nodes
1515  for (unsigned n=0;n<n_node;n++)
1516  {
1517  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(
1518  this->mesh_pt()->node_pt(n));
1519  for (unsigned i=0;i<n_dim;i++)
1520  {
1521  for(unsigned k=0;k<n_position_type;k++)
1522  {
1523  solid_nod_pt->pin_position(k,i);
1524  }
1525  }
1526  }
1527  }
1528  }
1529 
1530 
1531  ///\short Unpin all the field values and position unknowns (bit inefficient)
1532  void unpin_all()
1533  {
1534  // No need to do anything if there are no elements (in fact, we
1535  // probably never get here...)
1536  if (Problem::mesh_pt()->nelement()==0) return;
1537 
1538  //Loop over all the elements
1539  const unsigned n_element = Problem::mesh_pt()->nelement();
1540  for(unsigned e=0;e<n_element;++e)
1541  {
1542  //Cast the first element
1543  PROJECTABLE_ELEMENT * new_el_pt =
1544  dynamic_cast<PROJECTABLE_ELEMENT*>
1546  //Find the number of fields
1547  unsigned n_fields = new_el_pt->nfields_for_projection();
1548 
1549  //Now loop over all fields
1550  for(unsigned f=0;f<n_fields;f++)
1551  {
1552  //Extract the data and value for the field
1554  data=new_el_pt->data_values_of_field(f);
1555  //Now loop over all the data and unpin the appropriate values
1556  unsigned d_size=data.size();
1557  for(unsigned d=0;d<d_size;d++)
1558  {
1559  data[d].first->unpin(data[d].second);
1560  }
1561  }
1562  }
1563 
1564  /// Do we have a solid mesh?
1565  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>
1566  (Problem::mesh_pt()->element_pt(0));
1567  if (solid_el_pt!=0)
1568  {
1569  //Find number of nodes
1570  const unsigned n_node=this->mesh_pt()->nnode();
1571  //If no nodes then return
1572  if(n_node==0) {return;}
1573 
1574  //Read dimension and number of position values from the first node
1575  const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1576  const unsigned n_position_type =
1577  this->mesh_pt()->node_pt(0)->nposition_type();
1578 
1579  //Loop over the nodes
1580  for (unsigned n=0;n<n_node;n++)
1581  {
1582  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(
1583  this->mesh_pt()->node_pt(n));
1584  for (unsigned i=0;i<n_dim;i++)
1585  {
1586  for(unsigned k=0;k<n_position_type;k++)
1587  {
1588  solid_nod_pt->unpin_position(k,i);
1589  }
1590  }
1591  }
1592  }
1593  }
1594 
1595  /// \short Helper function to unpin the values of coordinate coord
1596  void unpin_dofs_of_coordinate(const unsigned &coord)
1597  {
1598  //Loop over the nodes
1599  const unsigned n_node = Problem::mesh_pt()->nnode();
1600  //If there are no nodes return immediately
1601  if(n_node==0) {return;}
1602 
1603  //Find the number of position values (should be the same for all nodes)
1604  const unsigned n_position_type =
1606 
1607  for(unsigned n=0;n<n_node;++n)
1608  {
1609  //Cache access to the Node as a solid node
1610  SolidNode* solid_nod_pt =
1611  static_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
1612  //Unpin all position types associated with the given coordinate
1613  for(unsigned k=0;k<n_position_type;++k)
1614  {
1615  solid_nod_pt->unpin_position(k,coord);
1616  }
1617  }
1618  }
1619 
1620  /// \short Helper function to unpin the values of coordinate coord
1621  void pin_dofs_of_coordinate(const unsigned &coord)
1622  {
1623  //Loop over the nodes
1624  const unsigned n_node = Problem::mesh_pt()->nnode();
1625  //If there are no nodes return immediately
1626  if(n_node==0) {return;}
1627 
1628  //Find the number of position values (should be the same for all nodes)
1629  const unsigned n_position_type =
1631 
1632  for(unsigned n=0;n<n_node;++n)
1633  {
1634  //Cache access to the Node as a solid node
1635  SolidNode* solid_nod_pt =
1636  static_cast<SolidNode*>(Problem::mesh_pt()->node_pt(n));
1637  //Unpin all position types associated with the given coordinate
1638  for(unsigned k=0;k<n_position_type;++k)
1639  {
1640  solid_nod_pt->pin_position(k,coord);
1641  }
1642  }
1643  }
1644 
1645 
1646  ///Helper function to unpin dofs of fld-th field
1647  void unpin_dofs_of_field(const unsigned &fld)
1648  {
1649  unsigned n_element = Problem::mesh_pt()->nelement();
1650  for(unsigned e=0;e<n_element;e++)
1651  {
1652  PROJECTABLE_ELEMENT * new_el_pt =
1653  dynamic_cast<PROJECTABLE_ELEMENT*>
1655 
1657  data=new_el_pt->data_values_of_field(fld);
1658 
1659  unsigned d_size=data.size();
1660  for(unsigned d=0;d<d_size;d++)
1661  {
1662  data[d].first->unpin(data[d].second);
1663  }
1664  }
1665  }
1666 
1667  ///Helper function to unpin dofs of fld-th field
1668  void pin_dofs_of_field(const unsigned &fld)
1669  {
1670  unsigned n_element = Problem::mesh_pt()->nelement();
1671  for(unsigned e=0;e<n_element;e++)
1672  {
1673  PROJECTABLE_ELEMENT * new_el_pt =
1674  dynamic_cast<PROJECTABLE_ELEMENT*>
1676 
1678  data=new_el_pt->data_values_of_field(fld);
1679 
1680  unsigned d_size=data.size();
1681  for(unsigned d=0;d<d_size;d++)
1682  {
1683  data[d].first->pin(data[d].second);
1684  }
1685  }
1686  }
1687 
1688  /// Helper function to set time level for projection
1689  void set_time_level_for_projection(const unsigned &time_level)
1690  {
1691  unsigned n_element = Problem::mesh_pt()->nelement();
1692  for(unsigned e=0;e<n_element;e++)
1693  {
1694  PROJECTABLE_ELEMENT * el_pt =
1695  dynamic_cast<PROJECTABLE_ELEMENT*>
1697 
1698  //Set what time we are dealing with
1699  el_pt->time_level_for_projection()=time_level;
1700  }
1701  }
1702 
1703  ///Set the coordinate for projection
1704  void set_coordinate_for_projection(const unsigned &i)
1705  {
1706  unsigned n_element = Problem::mesh_pt()->nelement();
1707  for(unsigned e=0;e<n_element;e++)
1708  {
1709  PROJECTABLE_ELEMENT * new_el_pt =
1710  dynamic_cast<PROJECTABLE_ELEMENT*>
1712 
1713  //Set that we are solving a projected coordinate problem
1714  new_el_pt->set_project_coordinates();
1715 
1716  //Set the projected coordinate
1717  new_el_pt->projected_coordinate()=i;
1718  }
1719  }
1720 
1721  ///Set the Lagrangian coordinate for projection
1723  {
1724  unsigned n_element = Problem::mesh_pt()->nelement();
1725  for(unsigned e=0;e<n_element;e++)
1726  {
1727  PROJECTABLE_ELEMENT * new_el_pt =
1728  dynamic_cast<PROJECTABLE_ELEMENT*>
1730 
1731  //Set that we are solving a projected Lagrangian coordinate problem
1732  new_el_pt->set_project_lagrangian();
1733 
1734  //Set the projected coordinate
1735  new_el_pt->projected_lagrangian_coordinate()=i;
1736  }
1737  }
1738 
1739  ///Set current field for projection
1740  void set_current_field_for_projection(const unsigned &fld)
1741  {
1742  unsigned n_element = Problem::mesh_pt()->nelement();
1743  for(unsigned e=0;e<n_element;e++)
1744  {
1745  PROJECTABLE_ELEMENT * new_el_pt =
1746  dynamic_cast<PROJECTABLE_ELEMENT*>
1748 
1749  //Set current field
1750  new_el_pt->projected_field()=fld;
1751  }
1752  }
1753 
1754  private:
1755 
1756  /// Backup for pin status of solid node's position Data
1758 
1759  /// Flag to suppress output during projection
1761 
1762  // Use an iterative solver for solving the system of equations
1764 
1765  // The iterative solver to solve the projection problem
1767 
1768  // The preconditioner for the solver
1770 
1771 };
1772 
1773 
1774 } //End of namespace
1775 
1776 #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'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:1704
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 newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8741
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
void unpin_dofs_of_coordinate(const unsigned &coord)
Helper function to unpin the values of coordinate coord.
Definition: projection.h:1596
double nodal_position_gen(const unsigned &n, const unsigned &k, const unsigned &i) const
Return the value of the k-th type of the i-th positional variable at the local node n...
Definition: elements.h:2237
void describe_local_dofs(std::ostream &out, const std::string &curr_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void setup_multi_domain_interaction(Problem *problem_pt, Mesh *const &mesh_pt, Mesh *const &external_mesh_pt, const unsigned &interaction_index=0)
Function to set up the one-way multi-domain interaction for problems where the meshes pointed to by m...
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
unsigned & time_level_for_projection()
Which history value are we projecting?
Definition: projection.h:627
IterativeLinearSolver * Iterative_solver_projection_pt
Definition: projection.h:1766
void unpin_position(const unsigned &i)
Unpin the nodal position.
Definition: nodes.h:1702
void pin_all()
Pin all the field values and position unknowns (bit inefficient)
Definition: projection.h:1465
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...
cstr elem_len * i
Definition: cfortran.h:607
void set_project_lagrangian()
Project (current and only values of) lagrangian coordinates.
Definition: projection.h:619
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
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:90
void pin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
Definition: projection.h:1668
Template-free Base class for projectable elements.
Definition: projection.h:59
The Problem class.
Definition: problem.h:152
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 enable_suppress_output_during_projection()
Suppress all output during projection phases.
Definition: projection.h:699
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Overloaded version of fill_in_contribution_to_residuals.
Definition: projection.h:187
A general Finite Element class.
Definition: elements.h:1271
unsigned & projected_field()
Field that is currently being projected.
Definition: projection.h:623
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void set_project_coordinates()
Project (history values of) coordintes.
Definition: projection.h:613
bool Shut_up_in_newton_solve
Boolean to indicate if all output is suppressed in Problem::newton_solve(). Defaults to false...
Definition: problem.h:2069
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void project(Mesh *base_mesh_pt, const bool &dont_project_positions=false)
Project from base into the problem's own mesh.
Definition: projection.h:724
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
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 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
void set_lagrangian_coordinate_for_projection(const unsigned &i)
Set the Lagrangian coordinate for projection.
Definition: projection.h:1722
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
void store_positions()
Helper function to store positions (the only things that have been set before doing projection...
Definition: projection.h:1387
OomphInfo oomph_info
An interface to the Trilinos AztecOO classes allowing it to be used as an Oomph-lib LinearSolver...
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting...
Definition: projection.h:76
static bool mpi_has_been_initialised()
return true if MPI has been initialised
e
Definition: cfortran.h:575
bool Output_during_projection_suppressed
Flag to suppress output during projection.
Definition: projection.h:1760
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Overloaded version of fill_in_contribution_to_jacobian.
Definition: projection.h:218
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void ignore_external_interaction_data()
Do not include any external interaction data when computing the element's Jacobian.
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
void unpin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
Definition: projection.h:1647
bool use_iterative_solver_for_projection()
Definition: projection.h:712
void disable_info_in_newton_solve()
Disable the output of information when in the newton solver.
Definition: problem.h:2812
Matrix-based diagonal preconditioner.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
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
bool add_external_geometric_data()
Are we including external geometric data in the element's Jacobian.
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1693
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
double & xi_gen(const unsigned &k, const unsigned &i)
Reference to the generalised Lagrangian position. `Type': k; 'Coordinate direction: i...
Definition: nodes.h:1759
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(...).
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
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 fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and the element Jacobian matrix (wrapper) ...
void enable_projection()
Backup the element'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:716
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
unsigned & projected_coordinate()
When projecting the history values of the nodal coordinates, this is the coordinate we'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:1621
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
void disable_doc_time()
Disable documentation of solve times.
unsigned & projected_lagrangian_coordinate()
When projecting the Lagrangian coordinates this is the coordinate that is being projected.
Definition: projection.h:637
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void unpin_all()
Unpin all the field values and position unknowns (bit inefficient)
Definition: projection.h:1532
double timer()
returns the time in seconds after some point in past
bool add_external_interaction_data()
Are we including external data in the element's Jacobian.
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...
friend class RefineableTetgenMesh
Definition: projection.h:687
void set_time_level_for_projection(const unsigned &time_level)
Helper function to set time level for projection.
Definition: projection.h:1689
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
void disable_use_iterative_solver_for_projection()
Disbales the use of an iterative solver for projection.
Definition: projection.h:720
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
void set_current_field_for_projection(const unsigned &fld)
Set current field for projection.
Definition: projection.h:1740
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
The conjugate gradient method.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1424
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
bool is_doc_time_enabled() const
Is documentation of solve times enabled?
void ignore_external_geometric_data()
Do not include any external geometric data when computing the element's Jacobian. This function shoul...
Vector< DenseMatrix< double > > Solid_backup
Backup for pin status of solid node's position Data.
Definition: projection.h:1757
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:705
void restore_positions()
Helper function to restore positions (the only things that have been set before doing projection...
Definition: projection.h:1428
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
unsigned ninteraction() const
Return the number of interactions in the element.
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1737
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
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 Problem_is_nonlinear
Boolean flag indicating if we're dealing with a linear or nonlinear Problem – if set to false the New...
Definition: problem.h:625
void enable_doc_time()
Enable documentation of solve times.
bool Use_iterative_solver_for_projection
Definition: projection.h:1763
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:3881
void include_external_geometric_data()
Do include external geometric data when computing the element's Jacobian. This function should be cal...
Projection_Type
Enumerated collection to specify which projection problem is to be solved.
Definition: projection.h:66
Unstructured refineable Triangle Mesh.
SolidFiniteElement class.
Definition: elements.h:3320
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:1769
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
void include_external_interaction_data()
Do include external geometric data when computing the element's Jacobian This function should be call...
A general mesh class.
Definition: mesh.h:74
void set_project_values()
Project (history values of) values.
Definition: projection.h:616
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.
Definition: projection.h:80