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-2022 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"
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
44
45// Use a preconditioner for the iterative solver
46#include "preconditioner.h"
48
49namespace 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
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.
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
187 {
188 // Do projection
189 if (Do_projection)
190 {
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 {
213 ELEMENT::describe_local_dofs(out, current_string);
214 }
215
216 /// Overloaded version of fill_in_contribution_to_jacobian
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
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 =
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
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;
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
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 {
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 {
648 }
649
650 /// When projecting the Lagrangian coordinates this is
651 /// the coordinate that is being projected
653 {
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>
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>
703 template<class FRIEND_PROJECTABLE_ELEMENT>
705 template<class FRIEND_PROJECTABLE_ELEMENT>
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?
758 {
759 // Create a Trilinos Solver
761 // Select CG as the linear solver
763 ->solver_type() = TrilinosAztecOOSolver::CG;
764 }
765 else
766 {
767 // Use CG to solve the projection problem
769 }
770
771 // Create the preconditioner
773 // Set the preconditioner for the solver
776
777 // Set CG as the linear solver
779
780#else
781 // Check whether the problem is distributed?
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
794
795 // Create the preconditioner
797 // Set the preconditioner for the solver
800
801 // Set CG as the linear solver
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 =
820 {
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
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();
907 PROJECTABLE_ELEMENT>(this, Problem::mesh_pt(), base_mesh_pt);
909 {
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
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 {
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
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
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*>(
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 {
1054 << "\n\n=============================================\n";
1055 oomph_info << "Projecting history values for coordinate " << i
1056 << std::endl;
1058 << "=============================================\n\n";
1059 }
1060
1061 // Set the coordinate for projection
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 {
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
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
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 {
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*>(
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 {
1176 << "\n\n=============================================\n";
1177 oomph_info << "Projecting history values for coordinate " << i
1178 << std::endl;
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 {
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
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*>(
1216
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
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*>(
1305
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 {
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 {
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)
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
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 =
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 =
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
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
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
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
1809
1810 // The preconditioner for the solver
1812 };
1813
1814
1815} // namespace oomph
1816
1817#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
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 set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
void include_external_geometric_data()
Do include external geometric data when computing the element's Jacobian. This function should be cal...
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 ...
bool add_external_geometric_data()
Are we including external geometric data in the element's Jacobian.
void ignore_external_interaction_data()
Do not include any external interaction data when computing the element's Jacobian.
bool add_external_interaction_data()
Are we including external data in the element's Jacobian.
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
void include_external_interaction_data()
Do include external geometric data when computing the element's Jacobian This function should be call...
void ignore_external_geometric_data()
Do not include any external geometric data when computing the element's Jacobian. This function shoul...
unsigned ninteraction() const
Return the number of interactions in the element.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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:4103
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
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:2349
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
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:227
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
bool is_doc_time_enabled() const
Is documentation of solve times enabled?
void disable_doc_time()
Disable documentation of solve times.
void enable_doc_time()
Enable documentation of solve times.
static bool mpi_has_been_initialised()
return true if MPI has been initialised
Matrix-based diagonal preconditioner.
A general mesh class.
Definition: mesh.h:67
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i). ‘Type’: k; Coordinate direction: i.
Definition: nodes.h:1126
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
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:1989
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
bool Shut_up_in_newton_solve
Boolean to indicate if all output is suppressed in Problem::newton_solve(). Defaults to false.
Definition: problem.h:2191
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
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:628
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
void disable_info_in_newton_solve()
Disable the output of information when in the newton solver.
Definition: problem.h:2941
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 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 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 ~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
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
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_field()
Field that is currently being projected.
Definition: projection.h:632
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
unsigned & projected_lagrangian_coordinate()
When projecting the Lagrangian coordinates this is the coordinate that is being projected.
Definition: projection.h:652
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
unsigned & projected_coordinate()
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting.
Definition: projection.h:645
Projection problem. This is created during the adaptation of unstructured meshes and it is assumed th...
Definition: projection.h:695
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
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
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
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
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:4137
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1877
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
void unpin_position(const unsigned &i)
Unpin the nodal position.
Definition: nodes.h:1829
double & xi_gen(const unsigned &k, const unsigned &i)
Reference to the generalised Lagrangian position. ‘Type’: k; 'Coordinate direction: i.
Definition: nodes.h:1902
An interface to the Trilinos AztecOO classes allowing it to be used as an Oomph-lib LinearSolver....
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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...
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...