fish_mesh.template.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_FISH_MESH_TEMPLATE_CC
27 #define OOMPH_FISH_MESH_TEMPLATE_CC
28 
29 #include "fish_mesh.template.h"
30 
31 
32 namespace oomph
33 {
34  //=================================================================
35  /// Constructor: Pass pointer to timestepper.
36  /// (defaults to (Steady) default timestepper defined in Mesh)
37  //=================================================================
38  template<class ELEMENT>
39  FishMesh<ELEMENT>::FishMesh(TimeStepper* time_stepper_pt)
40  {
41  // Mesh can only be built with 2D Qelements.
42  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
43 
44  // Fish back is a circle of radius 1, centred at (0.5,0.0)
45  double x_c = 0.5;
46  double y_c = 0.0;
47  double r_back = 1.0;
48  Back_pt = new Circle(x_c, y_c, r_back);
49 
50  // I've created the fishback -- I need to kill it.
51  Must_kill_fish_back = true;
52 
53  // Now build the mesh
54  build_mesh(time_stepper_pt);
55  }
56 
57 
58  //=================================================================
59  /// Constructor: Pass pointer GeomObject that defines
60  /// the fish's back and pointer to timestepper.
61  /// (defaults to (Steady) default timestepper defined in Mesh)
62  //=================================================================
63  template<class ELEMENT>
64  FishMesh<ELEMENT>::FishMesh(GeomObject* back_pt, TimeStepper* time_stepper_pt)
65  : Back_pt(back_pt)
66  {
67  // Mesh can only be built with 2D Qelements.
68  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
69 
70  // Back_pt has already been set....
71  Must_kill_fish_back = false;
72 
73  // Now build the mesh
74  build_mesh(time_stepper_pt);
75  }
76 
77 
78  //============================start_build_mesh=====================
79  /// Build the mesh, using the geometric object that
80  /// defines the fish's back.
81  //=================================================================
82  template<class ELEMENT>
83  void FishMesh<ELEMENT>::build_mesh(TimeStepper* time_stepper_pt)
84  {
85  // Mesh can only be built with 2D Qelements.
86  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
87 
88  // Build domain: Pass the pointer to the geometric object that
89  // represents the fish's back (pointed to by the FishMesh's
90  // private data member, Back_pt, and the values of the
91  // Lagrangian coordinate along this object at the "tail"
92  // and the "nose":
93  double xi_lo = 2.6;
94  double xi_hi = 0.4;
95  Domain_pt = new FishDomain(Back_pt, xi_lo, xi_hi);
96 
97  // Plot the domain? Keep this here in case we need to doc it
98  bool plot_it = false;
99  if (plot_it)
100  {
101  // Output the domain
102  std::ofstream some_file;
103 
104  // Number of plot points in each coordinate direction.
105  unsigned npts = 10;
106 
107  // Output domain
108  some_file.open("fish_domain.dat");
109  Domain_pt->output(some_file, npts);
110  Domain_pt->output_macro_element_boundaries(some_file, npts);
111  some_file.close();
112  }
113 
114  // Set the number of boundaries
115  set_nboundary(7);
116 
117 
118  // We will store boundary coordinates on the curvilinear boundaries
119  //(boundaries 0 and 4) along the fish's belly and its back.
120  Boundary_coordinate_exists[0] = true;
121  Boundary_coordinate_exists[4] = true;
122 
123  // Allocate the storage for the elements
124  unsigned nelem = 4;
125  Element_pt.resize(nelem);
126 
127  // Create dummy element so we can determine the number of nodes
128  ELEMENT* dummy_el_pt = new ELEMENT;
129 
130  // Read out the number of linear points in the element
131  unsigned n_node_1d = dummy_el_pt->nnode_1d();
132 
133  // Kill the element
134  delete dummy_el_pt;
135 
136  // Can now allocate the store for the nodes
137  unsigned nnodes_total =
138  (2 * (n_node_1d - 1) + 1) * (2 * (n_node_1d - 1) + 1);
139  Node_pt.resize(nnodes_total);
140 
141 
142  Vector<double> s(2), s_fraction(2);
143  Vector<double> r(2);
144 
145 
146  // Create elements and all nodes in element
147  //-----------------------------------------
148  // (ignore repetitions for now -- we'll clean them up later)
149  //----------------------------------------------------------
150  for (unsigned e = 0; e < nelem; e++)
151  {
152  // Create element
153  Element_pt[e] = new ELEMENT;
154 
155  // Loop over rows in y/s_1-direction
156  for (unsigned i1 = 0; i1 < n_node_1d; i1++)
157  {
158  // Loop over rows in x/s_0-direction
159  for (unsigned i0 = 0; i0 < n_node_1d; i0++)
160  {
161  // Local node number
162  unsigned j_local = i0 + i1 * n_node_1d;
163 
164  // Create the node and store pointer to it
165  Node* node_pt =
166  finite_element_pt(e)->construct_node(j_local, time_stepper_pt);
167 
168  // Work out the node's coordinates in the finite element's local
169  // coordinate system:
170  finite_element_pt(e)->local_fraction_of_node(j_local, s_fraction);
171 
172  s[0] = -1.0 + 2.0 * s_fraction[0];
173  s[1] = -1.0 + 2.0 * s_fraction[1];
174 
175  // Get the global position of the node from macro element mapping
176  Domain_pt->macro_element_pt(e)->macro_map(s, r);
177 
178  // Set the nodal position
179  node_pt->x(0) = r[0];
180  node_pt->x(1) = r[1];
181  }
182  }
183  } // end of loop over elements
184 
185 
186  // Kill repeated nodes and replace by pointers to nodes in appropriate
187  //---------------------------------------------------------------------
188  // neighbour. Also add node pointers to boundary arrays.
189  //------------------------------------------------------
190 
191  // Initialise number of global nodes
192  unsigned node_count = 0;
193 
194  // Check for error in node killing
195  bool stopit = false;
196 
197 
198  // Max tolerance for error in node killing
199  double Max_tol_in_node_killing = 1.0e-12;
200 
201  // First element: Lower body: all nodes survive
202  //---------------------------------------------
203  unsigned e = 0;
204 
205  // Loop over rows in y/s_1-direction
206  for (unsigned i1 = 0; i1 < n_node_1d; i1++)
207  {
208  // Loop over rows in x/s_0-direction
209  for (unsigned i0 = 0; i0 < n_node_1d; i0++)
210  {
211  // Local node number
212  unsigned j_local = i0 + i1 * n_node_1d;
213 
214  // No duplicate node: Copy new node across to mesh
215  Node_pt[node_count] = finite_element_pt(e)->node_pt(j_local);
216 
217  // Set boundaries:
218 
219  // If we're on the boundary need to convert the node
220  // into a boundary node
221  if ((i0 == 0) || (i1 == 0))
222  {
223  this->convert_to_boundary_node(Node_pt[node_count]);
224  }
225 
226  // Lower jaw: boundary 6
227  if (i0 == 0)
228  {
229  add_boundary_node(6, Node_pt[node_count]);
230  }
231 
232  // Belly: boundary 0
233  if (i1 == 0)
234  {
235  add_boundary_node(0, Node_pt[node_count]);
236 
237  // Set the boundary coordinate
238  Vector<double> zeta(1);
239  zeta[0] =
240  xi_lo + (xi_hi - xi_lo) * double(i0) / double(n_node_1d - 1);
241  Node_pt[node_count]->set_coordinates_on_boundary(0, zeta);
242  }
243 
244  // Increment node counter
245  node_count++;
246  }
247  }
248 
249 
250  // Lower fin: Western row of nodes is duplicate from element 0
251  //------------------------------------------------------------
252  e = 1;
253 
254  // Loop over rows in y/s_1-direction
255  for (unsigned i1 = 0; i1 < n_node_1d; i1++)
256  {
257  // Loop over rows in x/s_0-direction
258  for (unsigned i0 = 0; i0 < n_node_1d; i0++)
259  {
260  // Local node number
261  unsigned j_local = i0 + i1 * n_node_1d;
262 
263  // Has the node been killed?
264  bool killed = false;
265 
266  // First vertical row of nodes in s_1 direction get killed
267  // and re-directed to nodes in element 0
268  if (i0 == 0)
269  {
270  // Neighbour element
271  unsigned e_neigh = 0;
272 
273  // Node in neighbour element
274  unsigned i0_neigh = n_node_1d - 1;
275  unsigned i1_neigh = i1;
276  unsigned j_local_neigh = i0_neigh + i1_neigh * n_node_1d;
277 
278 
279  // Check:
280  for (unsigned i = 0; i < 2; i++)
281  {
282  double error = std::fabs(
283  finite_element_pt(e)->node_pt(j_local)->x(i) -
284  finite_element_pt(e_neigh)->node_pt(j_local_neigh)->x(i));
285  if (error > Max_tol_in_node_killing)
286  {
287  oomph_info << "Error in node killing for i " << i << " " << error
288  << std::endl;
289  stopit = true;
290  }
291  }
292 
293  // Kill node
294  delete finite_element_pt(e)->node_pt(j_local);
295  killed = true;
296 
297  // Set pointer to neighbour:
298  finite_element_pt(e)->node_pt(j_local) =
299  finite_element_pt(e_neigh)->node_pt(j_local_neigh);
300  }
301 
302 
303  // No duplicate node: Copy across to mesh
304  if (!killed)
305  {
306  // Copy the node across
307  Node_pt[node_count] = finite_element_pt(e)->node_pt(j_local);
308 
309  // If we're on a boundary turn the node into
310  // a boundary node
311  if ((i1 == 0) || (i0 == n_node_1d - 1))
312  {
313  this->convert_to_boundary_node(Node_pt[node_count]);
314  }
315 
316  // Increment node counter
317  node_count++;
318  }
319 
320  // Set boundaries:
321 
322  // Bottom of tail: boundary 1
323  if (i1 == 0)
324  {
325  add_boundary_node(1, finite_element_pt(e)->node_pt(j_local));
326  }
327 
328  // Vertical bit of tail: boundary 2
329  if (i0 == n_node_1d - 1)
330  {
331  add_boundary_node(2, finite_element_pt(e)->node_pt(j_local));
332  }
333  }
334  }
335 
336 
337  // Upper body: Southern row of nodes is duplicate from element 0
338  //--------------------------------------------------------------
339  e = 2;
340 
341  // Loop over rows in y/s_1-direction
342  for (unsigned i1 = 0; i1 < n_node_1d; i1++)
343  {
344  // Loop over rows in x/s_0-direction
345  for (unsigned i0 = 0; i0 < n_node_1d; i0++)
346  {
347  // Local node number
348  unsigned j_local = i0 + i1 * n_node_1d;
349 
350  // Has the node been killed?
351  bool killed = false;
352 
353  // First horizontal row of nodes in s_0 direction get killed
354  // and re-directed to nodes in element 0
355  if (i1 == 0)
356  {
357  // Neighbour element
358  unsigned e_neigh = 0;
359 
360  // Node in neighbour element
361  unsigned i0_neigh = i0;
362  unsigned i1_neigh = n_node_1d - 1;
363  unsigned j_local_neigh = i0_neigh + i1_neigh * n_node_1d;
364 
365  // Check:
366  for (unsigned i = 0; i < 2; i++)
367  {
368  double error = std::fabs(
369  finite_element_pt(e)->node_pt(j_local)->x(i) -
370  finite_element_pt(e_neigh)->node_pt(j_local_neigh)->x(i));
371  if (error > Max_tol_in_node_killing)
372  {
373  oomph_info << "Error in node killing for i " << i << " " << error
374  << std::endl;
375  stopit = true;
376  }
377  }
378 
379  // Kill node
380  delete finite_element_pt(e)->node_pt(j_local);
381  killed = true;
382 
383  // Set pointer to neighbour:
384  finite_element_pt(e)->node_pt(j_local) =
385  finite_element_pt(e_neigh)->node_pt(j_local_neigh);
386  }
387 
388  // No duplicate node: Copy across to mesh
389  if (!killed)
390  {
391  // Copy the old node across to the mesh
392  Node_pt[node_count] = finite_element_pt(e)->node_pt(j_local);
393 
394  // If we're on a boundary, convert the node into a boundary
395  // node. This will automatically update the entry in the mesh
396  if ((i1 == n_node_1d - 1) || (i0 == 0))
397  {
398  this->convert_to_boundary_node(Node_pt[node_count]);
399  }
400 
401  // Increment node counter
402  node_count++;
403  }
404 
405  // Set boundaries:
406 
407  // Back: boundary 4
408  if (i1 == n_node_1d - 1)
409  {
410  add_boundary_node(4, finite_element_pt(e)->node_pt(j_local));
411 
412  // Set the boundary coordinate
413  Vector<double> zeta(1);
414  zeta[0] =
415  xi_lo + (xi_hi - xi_lo) * double(i0) / double(n_node_1d - 1);
416  finite_element_pt(e)->node_pt(j_local)->set_coordinates_on_boundary(
417  4, zeta);
418  }
419 
420  // Upper jaw: boundary 5
421  if (i0 == 0)
422  {
423  add_boundary_node(5, finite_element_pt(e)->node_pt(j_local));
424  }
425  }
426  }
427 
428 
429  // Upper fin: Western/southern row of nodes is duplicate from element 2/1
430  //-----------------------------------------------------------------------
431  e = 3;
432 
433  // Loop over rows in y/s_1-direction
434  for (unsigned i1 = 0; i1 < n_node_1d; i1++)
435  {
436  // Loop over rows in x/s_0-direction
437  for (unsigned i0 = 0; i0 < n_node_1d; i0++)
438  {
439  // Local node number
440  unsigned j_local = i0 + i1 * n_node_1d;
441 
442  // Has the node been killed?
443  bool killed = false;
444 
445  // First vertical row of nodes in s_1 direction get killed
446  // and re-directed to nodes in element 2
447  if (i0 == 0)
448  {
449  // Neighbour element
450  unsigned e_neigh = 2;
451 
452  // Node in neighbour element
453  unsigned i0_neigh = n_node_1d - 1;
454  unsigned i1_neigh = i1;
455  unsigned j_local_neigh = i0_neigh + i1_neigh * n_node_1d;
456 
457 
458  // Check:
459  for (unsigned i = 0; i < 2; i++)
460  {
461  double error = std::fabs(
462  finite_element_pt(e)->node_pt(j_local)->x(i) -
463  finite_element_pt(e_neigh)->node_pt(j_local_neigh)->x(i));
464  if (error > Max_tol_in_node_killing)
465  {
466  oomph_info << "Error in node killing for i " << i << " " << error
467  << std::endl;
468  stopit = true;
469  }
470  }
471 
472  // Kill node
473  delete finite_element_pt(e)->node_pt(j_local);
474  killed = true;
475 
476  // Set pointer to neighbour:
477  finite_element_pt(e)->node_pt(j_local) =
478  finite_element_pt(e_neigh)->node_pt(j_local_neigh);
479  }
480 
481 
482  // First horizontal row of nodes in s_0 direction (apart from
483  // first node get killed and re-directed to nodes in element 1
484  if ((i0 != 0) && (i1 == 0))
485  {
486  // Neighbour element
487  unsigned e_neigh = 1;
488 
489  // Node in neighbour element
490  unsigned i0_neigh = i0;
491  unsigned i1_neigh = n_node_1d - 1;
492  unsigned j_local_neigh = i0_neigh + i1_neigh * n_node_1d;
493 
494 
495  // Check:
496  for (unsigned i = 0; i < 2; i++)
497  {
498  double error = std::fabs(
499  finite_element_pt(e)->node_pt(j_local)->x(i) -
500  finite_element_pt(e_neigh)->node_pt(j_local_neigh)->x(i));
501  if (error > Max_tol_in_node_killing)
502  {
503  oomph_info << "Error in node killing for i " << i << " " << error
504  << std::endl;
505  stopit = true;
506  }
507  }
508 
509  // Kill node
510  delete finite_element_pt(e)->node_pt(j_local);
511  killed = true;
512 
513  // Set pointer to neighbour:
514  finite_element_pt(e)->node_pt(j_local) =
515  finite_element_pt(e_neigh)->node_pt(j_local_neigh);
516  }
517 
518 
519  // No duplicate node: Copy across to mesh
520  if (!killed)
521  {
522  // Now copy the node across
523  Node_pt[node_count] = finite_element_pt(e)->node_pt(j_local);
524 
525  // If we're on the boundary, convert the node into
526  // a boundary node
527  if ((i1 == n_node_1d - 1) || (i0 == n_node_1d - 1))
528  {
529  this->convert_to_boundary_node(Node_pt[node_count]);
530  }
531 
532  // Increment node counter
533  node_count++;
534  }
535 
536  // Set boundaries:
537 
538  // To of tail: boundary 3
539  if (i1 == n_node_1d - 1)
540  {
541  add_boundary_node(3, finite_element_pt(e)->node_pt(j_local));
542  }
543 
544 
545  // Vertical bit of tail: boundary 2
546  if (i0 == n_node_1d - 1)
547  {
548  add_boundary_node(2, finite_element_pt(e)->node_pt(j_local));
549  }
550  }
551  }
552 
553  // Terminate if there's been an error
554  if (stopit)
555  {
556  std::ostringstream error_message;
557  error_message << "Error occured in node killing!\n";
558  error_message
559  << "Max. permitted difference in position of the two nodes\n "
560  << "that get 'merged' : " << Max_tol_in_node_killing << std::endl;
561 
562  throw OomphLibError(
563  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
564  }
565 
566 
567  // Loop over all elements and set macro element pointer
568  unsigned n_element = this->nelement();
569  for (unsigned e = 0; e < n_element; e++)
570  {
571  // Get pointer to element
572  FiniteElement* el_pt = this->finite_element_pt(e);
573 
574  // Set pointer to macro element to enable MacroElement-based
575  // remesh. Also enables the curvlinear boundaries
576  // of the mesh/domain get picked up during adaptive
577  // mesh refinement in derived classes.
578  el_pt->set_macro_elem_pt(this->Domain_pt->macro_element_pt(e));
579  }
580 
581 
582  // Setup boundary element lookup schemes
583  setup_boundary_element_info();
584 
585 
586  // Check the boundary coordinates
587 #ifdef PARANOID
588  {
589  Vector<double> zeta(1);
590  Vector<double> r(2);
591  bool stopit = false;
592  unsigned num_bound = nboundary();
593  for (unsigned ibound = 0; ibound < num_bound; ibound++)
594  {
595  if (boundary_coordinate_exists(ibound))
596  {
597  unsigned num_nod = nboundary_node(ibound);
598  for (unsigned inod = 0; inod < num_nod; inod++)
599  {
600  // Get the boundary coordinate
601  boundary_node_pt(ibound, inod)
602  ->get_coordinates_on_boundary(ibound, zeta);
603 
604  // Get position from wall object
605  Back_pt->position(zeta, r);
606 
607  // Flip it
608  if (ibound == 0) r[1] = -r[1];
609 
610  // Check:
611  for (unsigned i = 0; i < 2; i++)
612  {
613  double error =
614  std::fabs(r[i] - boundary_node_pt(ibound, inod)->x(i));
615  if (error > Max_tol_in_node_killing)
616  {
617  oomph_info << "Error in boundary coordinate for direction " << i
618  << " on boundary " << ibound << ":" << error
619  << std::endl;
620 
621  oomph_info << "x: " << r[0] << " "
622  << boundary_node_pt(ibound, inod)->x(0) << std::endl;
623 
624  oomph_info << "y: " << r[1] << " "
625  << boundary_node_pt(ibound, inod)->x(1) << std::endl
626  << std::endl;
627  stopit = true;
628  }
629  }
630  }
631  }
632  }
633 
634  // Terminate if there's been an error
635  if (stopit)
636  {
637  std::ostringstream error_message;
638  error_message << "Error occured in boundary coordinate setup!\n";
639  error_message << "Max. tolerance: " << Max_tol_in_node_killing
640  << std::endl;
641 
642  throw OomphLibError(error_message.str(),
643  OOMPH_CURRENT_FUNCTION,
644  OOMPH_EXCEPTION_LOCATION);
645  }
646  }
647 #endif
648  }
649 
650 
651  /// /////////////////////////////////////////////////////////////////////
652  /// /////////////////////////////////////////////////////////////////////
653  /// /////////////////////////////////////////////////////////////////////
654 
655 
656  //=========start_setup_adaptivity=========================================
657  /// Setup all the information that's required for spatial adaptivity:
658  /// Build quadtree forest.
659  //========================================================================
660  template<class ELEMENT>
662  {
663  // Setup quadtree forest
664  this->setup_quadtree_forest();
665 
666  } // end of setup_adaptivity
667 
668 
669  /// ////////////////////////////////////////////////////////////////////
670  /// ////////////////////////////////////////////////////////////////////
671  // AlgebraicElement fish-shaped mesh
672  /// ////////////////////////////////////////////////////////////////////
673  /// ////////////////////////////////////////////////////////////////////
674 
675 
676  //======================================================================
677  /// Setup algebraic update operation. Nodes are "suspended"
678  /// from the fish's back and the upper edge of the fin. Nodes
679  /// in the lower half are placed symmetrically.
680  //======================================================================
681  template<class ELEMENT>
683  {
684 #ifdef PARANOID
685  /// Pointer to algebraic element in lower body
686  AlgebraicElementBase* lower_body_pt =
687  dynamic_cast<AlgebraicElementBase*>(Mesh::element_pt(0));
688 
689  if (lower_body_pt == 0)
690  {
691  std::ostringstream error_message;
692  error_message << "Element in AlgebraicFishMesh must be\n"
693  << "derived from AlgebraicElementBase\n"
694  << "but it is of type: "
695  << typeid(Mesh::element_pt(0)).name() << std::endl;
696 
697  throw OomphLibError(
698  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
699  }
700 #endif
701 
702  // Read out the number of linear points in the element
703  unsigned n_p =
704  dynamic_cast<ELEMENT*>(FishMesh<ELEMENT>::Mesh::finite_element_pt(0))
705  ->nnode_1d();
706 
707  // Element 0: Lower body
708  //----------------------
709  {
710  unsigned ielem = 0;
711  FiniteElement* el_pt = Mesh::finite_element_pt(ielem);
712 
713  // Loop over rows in y/s_1-direction
714  for (unsigned i1 = 0; i1 < n_p; i1++)
715  {
716  // Loop over rows in x/s_0-direction
717  for (unsigned i0 = 0; i0 < n_p; i0++)
718  {
719  // Local node number
720  unsigned jnod = i0 + i1 * n_p;
721 
722  // One geometric object is involved in update operation
723  Vector<GeomObject*> geom_object_pt(1);
724  geom_object_pt[0] = this->Back_pt;
725 
726  // The update function requires three parameters:
727  Vector<double> ref_value(3);
728 
729  // First reference value: fractional x-position
730  ref_value[0] = double(i0) / double(n_p - 1);
731 
732  // Second reference value: fractional position along
733  // straight line from position on horizontal symmetry line to
734  // point on fish back
735  ref_value[1] = 1.0 - double(i1) / double(n_p - 1);
736 
737  // Third reference value: Sign (are we above or below the
738  // symmetry line?)
739  ref_value[2] = -1.0;
740 
741  // Setup algebraic update for node: Pass update information
742  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(jnod))
743  ->add_node_update_info(this->Lower_body, // enumerated ID
744  this, // mesh
745  geom_object_pt, // vector of geom objects
746  ref_value); // vector of ref. values
747  }
748  }
749  }
750 
751 
752  // Element 1: Lower fin
753  //---------------------
754  {
755  unsigned ielem = 1;
756  FiniteElement* el_pt = Mesh::finite_element_pt(ielem);
757 
758  // Loop over rows in y/s_1-direction
759  for (unsigned i1 = 0; i1 < n_p; i1++)
760  {
761  // Loop over rows in x/s_0-direction
762  for (unsigned i0 = 0; i0 < n_p; i0++)
763  {
764  // Local node number
765  unsigned jnod = i0 + i1 * n_p;
766 
767  // One geometric object is involved in update operation
768  Vector<GeomObject*> geom_object_pt(1);
769  geom_object_pt[0] = this->Back_pt;
770 
771  // The update function requires three parameters:
772  Vector<double> ref_value(3);
773 
774  // First reference value: fractional x-position
775  ref_value[0] = double(i0) / double(n_p - 1);
776 
777  // Second reference value: fractional position along
778  // straight line from position on horizontal symmetry line to
779  // point on fish back
780  ref_value[1] = 1.0 - double(i1) / double(n_p - 1);
781 
782  // Third reference value: Sign (are we above or below the
783  // symmetry line?)
784  ref_value[2] = -1.0;
785 
786  // Setup algebraic update for node: Pass update information
787  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(jnod))
788  ->add_node_update_info(this->Lower_fin, // enumerated ID
789  this, // mesh
790  geom_object_pt, // vector of geom objects
791  ref_value); // vector of ref. values
792  }
793  }
794  }
795 
796 
797  // Element 2: Upper body
798  //----------------------
799  {
800  unsigned ielem = 2;
801  FiniteElement* el_pt = Mesh::finite_element_pt(ielem);
802 
803  // Loop over rows in y/s_1-direction
804  for (unsigned i1 = 0; i1 < n_p; i1++)
805  {
806  // Loop over rows in x/s_0-direction
807  for (unsigned i0 = 0; i0 < n_p; i0++)
808  {
809  // Local node number
810  unsigned jnod = i0 + i1 * n_p;
811 
812  // One geometric object is involved in update operation
813  Vector<GeomObject*> geom_object_pt(1);
814  geom_object_pt[0] = this->Back_pt;
815 
816  // The update function requires three parameters:
817  Vector<double> ref_value(3);
818 
819  // First reference value: fractional x-position
820  ref_value[0] = double(i0) / double(n_p - 1);
821 
822  // Second reference value: fractional position along
823  // straight line from position on horizontal symmetry line to
824  // point on fish back
825  ref_value[1] = double(i1) / double(n_p - 1);
826 
827  // Third reference value: Sign (are we above or below the
828  // symmetry line?)
829  ref_value[2] = 1.0;
830 
831  // Setup algebraic update for node: Pass update information
832  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(jnod))
833  ->add_node_update_info(this->Upper_body, // enumerated ID
834  this, // mesh
835  geom_object_pt, // vector of geom objects
836  ref_value); // vector of ref. values
837  }
838  }
839  }
840 
841 
842  // Element 3: Upper fin
843  //---------------------
844  {
845  unsigned ielem = 3;
846  FiniteElement* el_pt = Mesh::finite_element_pt(ielem);
847 
848  // Loop over rows in y/s_1-direction
849  for (unsigned i1 = 0; i1 < n_p; i1++)
850  {
851  // Loop over rows in x/s_0-direction
852  for (unsigned i0 = 0; i0 < n_p; i0++)
853  {
854  // Local node number
855  unsigned jnod = i0 + i1 * n_p;
856 
857  // One geometric object is involved in update operation
858  Vector<GeomObject*> geom_object_pt(1);
859  geom_object_pt[0] = this->Back_pt;
860 
861  // The update function requires three parameters:
862  Vector<double> ref_value(3);
863 
864  // First reference value: fractional x-position
865  ref_value[0] = double(i0) / double(n_p - 1);
866 
867  // Second reference value: fractional position along
868  // straight line from position on horizontal symmetry line to
869  // point on fish back
870  ref_value[1] = double(i1) / double(n_p - 1);
871 
872  // Third reference value: Sign (are we above or below the
873  // symmetry line?)
874  ref_value[2] = 1.0;
875 
876  // Setup algebraic update for node: Pass update information
877  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(jnod))
878  ->add_node_update_info(this->Upper_fin, // enumerated ID
879  this, // mesh
880  geom_object_pt, // vector of geom objects
881  ref_value); // vector of ref. values
882  }
883  }
884  }
885  }
886 
887 
888  //======================================================================
889  /// Algebraic update function: Update in (upper or lower) body
890  /// according to wall shape at time level t (t=0: present; t>0: previous)
891  //======================================================================
892  template<class ELEMENT>
894  AlgebraicNode*& node_pt)
895  {
896  // Pointer to geometric object that represents the fish back:
897  GeomObject* back_pt = node_pt->geom_object_pt(unsigned(0));
898 
899  // Fixed reference value: x-position of mouth
900  double x_mouth = this->Domain_pt->x_mouth();
901 
902  // Fixed reference value: Lagrangian coordinate of point
903  // over mouth
904  double zeta_mouth = this->Domain_pt->xi_nose();
905 
906  // Fixed reference value: Lagrangian coordinate of point
907  // near tail
908  double zeta_near_tail = this->Domain_pt->xi_tail();
909 
910  // First reference value: fractional x-position
911  double fract_x = node_pt->ref_value(unsigned(0));
912 
913  // Second reference value: fractional position along
914  // straight line from position on horizontal symmetry line to
915  // point on fish back
916  double fract_y = node_pt->ref_value(1);
917 
918  // Third reference value: Sign (are we above or below the
919  // symmetry line?)
920  double sign = node_pt->ref_value(2);
921 
922  // Get position on fish back
923  Vector<double> zeta(back_pt->nlagrangian());
924  zeta[0] = zeta_mouth + fract_x * (zeta_near_tail - zeta_mouth);
925  Vector<double> r_back(back_pt->ndim());
926  back_pt->position(t, zeta, r_back);
927 
928  // Get position of point on fish back near tail
929  zeta[0] = zeta_near_tail;
930  Vector<double> r_near_tail(back_pt->ndim());
931  back_pt->position(t, zeta, r_near_tail);
932 
933  // Get position on symmetry line
934  Vector<double> r_sym(2);
935  r_sym[0] = x_mouth + fract_x * (r_near_tail[0] - x_mouth);
936  r_sym[1] = 0.0;
937 
938  // Assign new nodal coordinate
939  node_pt->x(t, 0) = r_sym[0] + fract_y * (r_back[0] - r_sym[0]);
940  node_pt->x(t, 1) = sign * (r_sym[1] + fract_y * (r_back[1] - r_sym[1]));
941  }
942 
943 
944  //======================================================================
945  /// Algebraic update function: Update in (upper or lower) fin
946  /// according to wall shape at time level t (t=0: present; t>0: previous)
947  //======================================================================
948  template<class ELEMENT>
950  AlgebraicNode*& node_pt)
951  {
952  // Pointer to geometric object that represents the fish back:
953  GeomObject* back_pt = node_pt->geom_object_pt(unsigned(0));
954 
955  // Fixed reference value: End coordinate on fish back
956  double zeta_wall = this->Domain_pt->xi_tail();
957 
958  // Fixed reference value: x-position of end of tail
959  double x_tail = this->Domain_pt->x_fin();
960 
961  // Fixed reference value: y-position of end of tail
962  double y_tail = this->Domain_pt->y_fin();
963 
964  // First reference value: fractional position in x-direction
965  double fract_x = node_pt->ref_value(unsigned(0));
966 
967  // Second reference value: fractional position along
968  // vertical line from position on horizontal symmetry line to
969  // point on upper end of tail
970  double fract_y = node_pt->ref_value(1);
971 
972  // Third reference value: Sign (are we above or below the
973  // symmetry line?)
974  double sign = node_pt->ref_value(2);
975 
976  // Get position on fish back
977  Vector<double> zeta(back_pt->nlagrangian());
978  zeta[0] = zeta_wall;
979  Vector<double> r_back(back_pt->ndim());
980  back_pt->position(t, zeta, r_back);
981 
982  // y-position on top edge of fin:
983  double y_fin_edge = r_back[1] + fract_x * (y_tail - r_back[1]);
984 
985  // Assign new nodal coordinate
986  node_pt->x(t, 0) = r_back[0] + fract_x * (x_tail - r_back[0]);
987  node_pt->x(t, 1) = sign * fract_y * y_fin_edge;
988  }
989 
990 } // namespace oomph
991 #endif
void setup_algebraic_node_update()
Setup algebraic update operation for all nodes (separate function because this task needs to be perfo...
void node_update_in_body(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for nodes in upper/lower body.
void node_update_in_fin(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for nodes in upper/lower fin.
Fish shaped domain, represented by four MacroElements. Shape is parametrised by GeomObject that repre...
Definition: fish_domain.h:43
Fish shaped mesh. The geometry is defined by the Domain object FishDomain.
FishMesh(TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to timestepper (defaults to the (Steady) default timestepper defined in Mes...
void build_mesh(TimeStepper *time_stepper_pt)
Build the mesh, using the geometric object identified by Back_pt.
bool Must_kill_fish_back
Do I need to kill the fish back geom object?
void setup_adaptivity()
Setup all the information that's required for spatial adaptivity: Set pointers to macro elements and ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...