collapsible_channel_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_COLLAPSIBLE_CHANNEL_MESH_TEMPLATE_CC
31 #define OOMPH_COLLAPSIBLE_CHANNEL_MESH_TEMPLATE_CC
32 
33 //Include the headers file for collapsible channel
35 
36 
37 namespace oomph
38 {
39 
40 //========================================================================
41 /// Constructor: Pass number of elements in upstream/collapsible/downstream
42 /// segment and across the channel; lengths of upstream/collapsible/downstream
43 /// segments and width of channel, pointer to GeomObject that defines
44 /// the collapsible segment and pointer to TimeStepper (defaults to the
45 /// default timestepper, Steady).
46 //========================================================================
47 template <class ELEMENT>
49  const unsigned& nup,
50  const unsigned& ncollapsible,
51  const unsigned& ndown,
52  const unsigned& ny,
53  const double& lup,
54  const double& lcollapsible,
55  const double& ldown,
56  const double& ly,
57  GeomObject* wall_pt,
58  TimeStepper* time_stepper_pt)
59  : SimpleRectangularQuadMesh<ELEMENT>(nup+ncollapsible+ndown, ny,
60  lup+lcollapsible+ldown, ly,
61  time_stepper_pt),
62  Nup(nup), Ncollapsible(ncollapsible), Ndown(ndown), Ny(ny),
63  Wall_pt(wall_pt)
64 {
65  // Mesh can only be built with 2D Qelements.
66  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
67 
68  //Create the CollapsibleChannelDomain with the wall represented
69  //by the geometric object
70  Domain_pt = new CollapsibleChannelDomain(nup, ncollapsible, ndown, ny,
71  lup,lcollapsible, ldown, ly,
72  wall_pt);
73 
74  // Total number of (macro/finite)elements
75  unsigned nmacro=(nup+ncollapsible+ndown)*ny;
76 
77  // Loop over all elements and set macro element pointer
78  for (unsigned e=0;e<nmacro;e++)
79  {
80  // Get pointer to finite element
81  FiniteElement* el_pt=this->finite_element_pt(e);
82 
83  // Set pointer to macro element so the curvlinear boundaries
84  // of the mesh/domain can get picked up during adaptive
85  // mesh refinement
86  el_pt->set_macro_elem_pt(this->Domain_pt->macro_element_pt(e));
87  }
88 
89  // Update the nodal positions corresponding to their
90  // macro element representations.
91  this->node_update();
92 
93  // Update the boundary numbering scheme and set boundary coordinate
94  //-----------------------------------------------------------------
95 
96  // (Note: The original SimpleRectangularQuadMesh had four boundaries.
97  // We need to overwrite the boundary lookup scheme for the current
98  // mesh so that the collapsible segment becomes identifiable).
99  // While we're doing this, we're also setting up a boundary
100  // coordinate for the nodes located on the collapsible segment.
101  // The boundary coordinate can be used to setup FSI.
102 
103  // How many boundaries does the mesh have now?
104  unsigned nbound=this->nboundary();
105  for (unsigned b=0;b<nbound;b++)
106  {
107  // Remove all nodes on this boundary from the mesh's lookup scheme
108  // and also delete the reverse lookup scheme held by the nodes
109  this->remove_boundary_nodes(b);
110  }
111 
112 #ifdef PARANOID
113  // Sanity check
114  unsigned nnod=this->nnode();
115  for (unsigned j=0;j<nnod;j++)
116  {
117  if (this->node_pt(j)->is_on_boundary())
118  {
119  std::ostringstream error_message;
120  error_message << "Node " << j << "is still on boundary " << std::endl;
121 
122  throw OomphLibError(error_message.str(),
123  OOMPH_CURRENT_FUNCTION,
124  OOMPH_EXCEPTION_LOCATION);
125  }
126  }
127 #endif
128 
129  //Change the numbers of boundaries
130  this->set_nboundary(6);
131 
132  // Get the number of nodes along the element edge from first element
133  unsigned nnode_1d=this->finite_element_pt(0)->nnode_1d();
134 
135  // Vector of Lagrangian coordinates used as boundary coordinate
136  Vector<double> zeta(1);
137 
138  // Index of first element underneath the collapsible bit
139  unsigned first_collapsible=(ny-1)*(nup+ncollapsible+ndown)+nup;
140 
141  // Zeta increment over elements (used for assignment of
142  // boundary coordinate)
143  double dzeta= lcollapsible/double(ncollapsible);
144 
145  // Manually loop over the elements near the boundaries and
146  // assign nodes to boundaries. Also set up boundary coordinate
147  unsigned nelem=this->nelement();
148  for (unsigned e=0;e<nelem;e++)
149  {
150  // Bottom row of elements
151  if (e<nup+ncollapsible+ndown)
152  {
153  for (unsigned i=0;i<nnode_1d;i++)
154  {
155  this->add_boundary_node(0, this->finite_element_pt(e)->node_pt(i));
156  }
157  }
158  // Upstream upper rigid bit
159  if ((e>(ny-1)*(nup+ncollapsible+ndown)-1)&&
160  (e<(ny-1)*(nup+ncollapsible+ndown)+nup))
161  {
162  for (unsigned i=0;i<nnode_1d;i++)
163  {
164  this->add_boundary_node(4,
165  this->finite_element_pt(e)->node_pt((nnode_1d-1)*nnode_1d+i));
166  }
167  }
168  // Collapsible bit
169  if ((e>(ny-1)*(nup+ncollapsible+ndown)+nup-1)&&
170  (e<(ny-1)*(nup+ncollapsible+ndown)+nup+ncollapsible))
171  {
172  for (unsigned i=0;i<nnode_1d;i++)
173  {
174  this->add_boundary_node(3,
175  this->finite_element_pt(e)->node_pt((nnode_1d-1)*nnode_1d+i));
176 
177  // What column of elements are we in?
178  unsigned ix=e-first_collapsible;
179 
180  // Zeta coordinate
181  zeta[0]=double(ix)*dzeta+double(i)*dzeta/double(nnode_1d-1);
182 
183  // Set boundary coordinate
184  this->finite_element_pt(e)->node_pt((nnode_1d-1)*nnode_1d+i)->
185  set_coordinates_on_boundary(3,zeta);
186  }
187  }
188  // Downstream upper rigid bit
189  if ((e>(ny-1)*(nup+ncollapsible+ndown)+nup+ncollapsible-1)&&
190  (e<ny*(nup+ncollapsible+ndown)))
191  {
192  for (unsigned i=0;i<nnode_1d;i++)
193  {
194  this->add_boundary_node(2,
195  this->finite_element_pt(e)->node_pt((nnode_1d-1)*nnode_1d+i));
196  }
197  }
198  // Left end
199  if (e%(nup+ncollapsible+ndown)==0)
200  {
201  for (unsigned i=0;i<nnode_1d;i++)
202  {
203  this->add_boundary_node(5,
204  this->finite_element_pt(e)->node_pt(i*nnode_1d));
205  }
206  }
207  // Right end
208  if (e%(nup+ncollapsible+ndown)==(nup+ncollapsible+ndown)-1)
209  {
210  for (unsigned i=0;i<nnode_1d;i++)
211  {
212  this->add_boundary_node(1,
213  this->finite_element_pt(e)->node_pt((i+1)*nnode_1d-1));
214  }
215  }
216  }
217 
218  // Re-setup lookup scheme that establishes which elements are located
219  // on the mesh boundaries (doesn't need to be wiped)
220  this->setup_boundary_element_info();
221 
222  //We have only bothered to parametrise boundary 3
223  this->Boundary_coordinate_exists[3] = true;
224 }
225 
226 
227 
228 
229 ///////////////////////////////////////////////////////////////////////////
230 ///////////////////////////////////////////////////////////////////////////
231 ///////////////////////////////////////////////////////////////////////////
232 
233 
234 
235 //=================================================================
236 /// Perform algebraic mesh update at time level t (t=0: present;
237 /// t>0: previous)
238 //=================================================================
239 template<class ELEMENT>
241  const unsigned& t, AlgebraicNode*& node_pt)
242 {
243 
244 
245 #ifdef PARANOID
246  // We're updating the nodal positions (!) at time level t
247  // and determine them by evaluating the wall GeomObject's
248  // position at that gime level. I believe this only makes sense
249  // if the t-th history value in the positional timestepper
250  // actually represents previous values (rather than some
251  // generalised quantity). Hence if this function is called with
252  // t>nprev_values(), we issue a warning and terminate the execution.
253  // It *might* be possible that the code still works correctly
254  // even if this condition is violated (e.g. if the GeomObject's
255  // position() function returns the appropriate "generalised"
256  // position value that is required by the timestepping scheme but it's
257  // probably worth flagging this up and forcing the user to manually switch
258  // off this warning if he/she is 100% sure that this is kosher.
259  if (t>node_pt->position_time_stepper_pt()->nprev_values())
260  {
261  std::string error_message =
262  "Trying to update the nodal position at a time level";
263  error_message += "beyond the number of previous values in the nodes'";
264  error_message += "position timestepper. This seems highly suspect!";
265  error_message += "If you're sure the code behaves correctly";
266  error_message += "in your application, remove this warning ";
267  error_message += "or recompile with PARNOID switched off.";
268 
269  std::string function_name =
270  "AlgebraicCollapsibleChannelMesh::";
271  function_name += "algebraic_node_update()";
272 
273  throw OomphLibError(error_message,
274  OOMPH_CURRENT_FUNCTION,
275  OOMPH_EXCEPTION_LOCATION);
276  }
277 #endif
278 
279  // Extract references for update by copy construction
280  Vector<double> ref_value(node_pt->vector_ref_value());
281 
282  // First reference value: Original x-position. Used as the start point
283  // for the lines connecting the nodes in the vertical direction
284  double x_bottom=ref_value[0];
285 
286  // Second reference value: Fractional position along
287  // straight line from the bottom (at the original x position)
288  // to the reference point on the upper wall
289  double fract=ref_value[1];
290 
291  // Third reference value: Reference local coordinate
292  // in GeomObject that represents the upper wall (local coordinate
293  // in finite element if the wall GeomObject is a finite element mesh)
294  Vector<double> s(1);
295  s[0]=ref_value[2];
296 
297  // Fourth reference value: zeta coordinate on the upper wall
298  // If the wall is a simple GeomObject, zeta[0]=s[0]
299  // but if it's a compound GeomObject (e.g. a finite element mesh)
300  // zeta scales during mesh refinement, whereas s[0] and the
301  // pointer to the geom object have to be re-computed.
302  // double zeta=ref_value[3]; // not needed here
303 
304  // Extract geometric objects for update by copy construction
305  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
306 
307  // Pointer to actual wall geom object (either the same as the wall object
308  // or the pointer to the actual finite element)
309  GeomObject* geom_obj_pt=geom_object_pt[0];
310 
311  // Get position vector to wall at previous timestep t!
312  Vector<double> r_wall(2);
313  geom_obj_pt->position(t,s,r_wall);
314 
315  // Assign new nodal coordinate
316  node_pt->x(t,0)=x_bottom+fract*(r_wall[0]-x_bottom);
317  node_pt->x(t,1)=fract*r_wall[1];
318 
319 }
320 
321 
322 
323 
324 
325 //=====start_setup=================================================
326 /// Setup algebraic mesh update -- assumes that mesh has
327 /// initially been set up with a flush upper wall
328 //=================================================================
329 template<class ELEMENT>
331 {
332  // Shorthand for some geometric data:
333  double l_up=this->domain_pt()->l_up();
334  double l_collapsible=this->domain_pt()->l_collapsible();
335 
336  // Loop over all nodes in mesh
337  unsigned nnod=this->nnode();
338  for (unsigned j=0;j<nnod;j++)
339  {
340  // Get pointer to node -- recall that that Mesh::node_pt(...) has been
341  // overloaded in the AlgebraicMesh class to return a pointer to
342  // an AlgebraicNode.
343  AlgebraicNode* nod_pt=node_pt(j);
344 
345  // Get coordinates
346  double x=nod_pt->x(0);
347  double y=nod_pt->x(1);
348 
349  // Check if it's in the collapsible part:
350  if ( (x>=l_up) && (x<=(l_up+l_collapsible)) )
351  {
352 
353  // Get zeta coordinate on the undeformed wall
354  Vector<double> zeta(1);
355  zeta[0]=x-l_up;
356 
357  // Get pointer to geometric (sub-)object and Lagrangian coordinate
358  // on that sub-object. For a wall that is represented by
359  // a single geom object, this simply returns the input.
360  // If the geom object consists of sub-objects (e.g.
361  // if it is a finite element mesh representing a wall,
362  // then we'll obtain the pointer to the finite element
363  // (in its incarnation as a GeomObject) and the
364  // local coordinate in that element.
365  GeomObject* geom_obj_pt;
366  Vector<double> s(1);
367  this->Wall_pt->locate_zeta(zeta,geom_obj_pt,s);
368 
369  // Get position vector to wall:
370  Vector<double> r_wall(2);
371  geom_obj_pt->position(s,r_wall);
372 
373  // Sanity check: Confirm that the wall is in its undeformed position
374 #ifdef PARANOID
375  if ((std::fabs(r_wall[0]-x)>1.0e-15)&&(std::fabs(r_wall[1]-y)>1.0e-15))
376  {
377  std::ostringstream error_stream;
378  error_stream
379  << "Wall must be in its undeformed position when\n"
380  << "algebraic node update information is set up!\n "
381  << "x-discrepancy: " << std::fabs(r_wall[0]-x) << std::endl
382  << "y-discrepancy: " << std::fabs(r_wall[1]-y) << std::endl;
383 
384  throw OomphLibError(
385  error_stream.str(),
386  OOMPH_CURRENT_FUNCTION,
387  OOMPH_EXCEPTION_LOCATION);
388  }
389 #endif
390 
391 
392 
393  // One geometric object is involved in update operation
394  Vector<GeomObject*> geom_object_pt(1);
395 
396  // The actual geometric object (If the wall is simple GeomObject
397  // this is the same as Wall_pt; if it's a compound GeomObject
398  // this points to the sub-object)
399  geom_object_pt[0]=geom_obj_pt;
400 
401  // The update function requires four parameters:
402  Vector<double> ref_value(4);
403 
404  // First reference value: Original x-position
405  ref_value[0]=r_wall[0];
406 
407  // Second reference value: fractional position along
408  // straight line from the bottom (at the original x position)
409  // to the point on the wall)
410  ref_value[1]=y/r_wall[1];
411 
412  // Third reference value: Reference local coordinate
413  // in wall element (local coordinate in FE if we're dealing
414  // with a wall mesh)
415  ref_value[2]=s[0];
416 
417  // Fourth reference value: zeta coordinate on wall
418  // If the wall is a simple GeomObject, zeta[0]=s[0]
419  // but if it's a compound GeomObject (e.g. a finite element mesh)
420  // zeta scales during mesh refinement, whereas s[0] and the
421  // pointer to the geom object have to be re-computed.
422  ref_value[3]=zeta[0];
423 
424  // Setup algebraic update for node: Pass update information
425  nod_pt->add_node_update_info(
426  this, // mesh
427  geom_object_pt, // vector of geom objects
428  ref_value); // vector of ref. values
429  }
430 
431  }
432 
433 } //end of setup_algebraic_node_update
434 
435 
436 
437 
438 ////////////////////////////////////////////////////////////////////
439 ////////////////////////////////////////////////////////////////////
440 ////////////////////////////////////////////////////////////////////
441 
442 
443 
444 
445 
446 
447 //========start_update_node_update=================================
448 /// Update the geometric references that are used
449 /// to update node after mesh adaptation.
450 //=================================================================
451 template<class ELEMENT>
453  AlgebraicNode*& node_pt)
454 {
455  // Extract reference values for node update by copy construction
456  Vector<double> ref_value(node_pt->vector_ref_value());
457 
458  // First reference value: Original x-position
459  // double x_bottom=ref_value[0]; // not needed here
460 
461  // Second reference value: fractional position along
462  // straight line from the bottom (at the original x position)
463  // to the point on the wall)
464  // double fract=ref_value[1]; // not needed here
465 
466  // Third reference value: Reference local coordinate
467  // in GeomObject (local coordinate in finite element if the wall
468  // GeomObject is a finite element mesh)
469  // Vector<double> s(1);
470  // s[0]=ref_value[2]; // This needs to be re-computed!
471 
472  // Fourth reference value: intrinsic coordinate on the (possibly
473  // compound) wall.
474  double zeta=ref_value[3];
475 
476  // Extract geometric objects for update by copy construction
477  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
478 
479  // Pointer to actual wall geom object (either the same as wall object
480  // or the pointer to the actual finite element)
481  //GeomObject* geom_obj_pt=geom_object_pt[0]; // This needs to be re-computed!
482 
483  // Get zeta coordinate on wall (as vector)
484  Vector<double> zeta_wall(1);
485  zeta_wall[0]=zeta;
486 
487  // Get pointer to geometric (sub-)object and Lagrangian coordinate
488  // on that sub-object. For a wall that is represented by
489  // a single geom object, this simply returns the input.
490  // If the geom object consists of sub-objects (e.g.
491  // if it is a finite element mesh representing a wall,
492  // then we'll obtain the pointer to the finite element
493  // (in its incarnation as a GeomObject) and the
494  // local coordinate in that element.
495  Vector<double> s(1);
496  GeomObject* geom_obj_pt;
497  this->Wall_pt->locate_zeta(zeta_wall,geom_obj_pt,s);
498 
499  // Update the pointer to the (sub-)GeomObject within which the
500  // reference point is located. (If the wall is simple GeomObject
501  // this is the same as Wall_pt; if it's a compound GeomObject
502  // this points to the sub-object)
503  geom_object_pt[0]=geom_obj_pt;
504 
505  // First reference value: Original x-position
506  // ref_value[0]=r_wall[0]; // unchanged
507 
508  // Second reference value: fractional position along
509  // straight line from the bottom (at the original x position)
510  // to the point on the wall)
511  // ref_value[1]=y/r_wall[1]; // unchanged
512 
513  // Update third reference value: Reference local coordinate
514  // in wall element (local coordinate in FE if we're dealing
515  // with a wall mesh)
516  ref_value[2]=s[0];
517 
518  // Fourth reference value: zeta coordinate on wall
519  // If the wall is a simple GeomObject, zeta[0]=s[0]
520  // but if it's a compound GeomObject (e.g. a finite element mesh)
521  // zeta scales during mesh refinement, whereas s[0] and the
522  // pointer to the geom object have to be re-computed.
523  // ref_value[3]=zeta[0]; //unchanged
524 
525 
526  // Kill the existing node update info
527  node_pt->kill_node_update_info();
528 
529  // Setup algebraic update for node: Pass update information
530  node_pt->add_node_update_info(
531  this, // mesh
532  geom_object_pt, // vector of geom objects
533  ref_value); // vector of ref. values
534 
535 }
536 
537 
538 
539 }
540 #endif
CollapsibleChannelMesh(const unsigned &nup, const unsigned &ncollapsible, const unsigned &ndown, const unsigned &ny, const double &lup, const double &lcollapsible, const double &ldown, const double &ly, GeomObject *wall_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in upstream/collapsible/ downstream segment and across the chann...
CollapsibleChannelDomain * Domain_pt
Pointer to domain.
void update_node_update(AlgebraicNode *&node_pt)
Update the node update data for specified node following any mesh adapation.
void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)
Update nodal position at time level t (t=0: present; t>0: previous)
void setup_algebraic_node_update()
Function to setup the algebraic node update.