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