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-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#ifndef OOMPH_FISH_MESH_TEMPLATE_CC
27#define OOMPH_FISH_MESH_TEMPLATE_CC
28
29#include "fish_mesh.template.h"
30
31
32namespace oomph
33{
34 //=================================================================
35 /// Constructor: Pass pointer to timestepper.
36 /// (defaults to (Steady) default timestepper defined in Mesh)
37 //=================================================================
38 template<class ELEMENT>
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>
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>
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;
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;
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;
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;
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
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
////////////////////////////////////////////////////////////////////
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.
////////////////////////////////////////////////////////////////////
GeomObject * geom_object_pt(const unsigned &i)
Return pointer to i-th geometric object involved in default (usually first) update function.
double ref_value(const unsigned &i)
Return i-th reference value involved in default (usually first) update function.
////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:873
A general Finite Element class.
Definition: elements.h:1313
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:1872
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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?
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:177
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:171
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
void setup_adaptivity()
Setup all the information that's required for spatial adaptivity: Set pointers to macro elements and ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...