unstructured_two_d_fluid.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 // Driver for flow past odd-shaped obstacle -- domain meshed with triangle.
31 // This is a warm-up problem for an unstructured fsi problem.
32 
33 //Generic includes
34 #include "generic.h"
35 #include "navier_stokes.h"
36 #include "solid.h"
37 #include "constitutive.h"
38 
39 // The mesh
40 #include "meshes/triangle_mesh.h"
41 
42 
43 using namespace std;
44 
45 using namespace oomph;
46 
47 
48 //===========start_mesh====================================================
49 /// Triangle-based mesh upgraded to become a (pseudo-) solid mesh
50 //=========================================================================
51 template<class ELEMENT>
52 class ElasticTriangleMesh : public virtual TriangleMesh<ELEMENT>,
53  public virtual SolidMesh
54 {
55 
56 public:
57 
58  /// \short Constructor:
59  ElasticTriangleMesh(const std::string& node_file_name,
60  const std::string& element_file_name,
61  const std::string& poly_file_name,
62  TimeStepper* time_stepper_pt=
63  &Mesh::Default_TimeStepper) :
64  TriangleMesh<ELEMENT>(node_file_name,element_file_name,
65  poly_file_name, time_stepper_pt)
66  {
67  //Assign the Lagrangian coordinates
68  set_lagrangian_nodal_coordinates();
69 
70  //Identify the domain boundaries
71  this->identify_boundaries();
72  }
73 
74 
75  ///\short Function used to identify the domain boundaries
77  {
78  // We will have three boundaries
79  this->set_nboundary(3);
80 
81  unsigned n_node=this->nnode();
82  for (unsigned j=0;j<n_node;j++)
83  {
84  Node* nod_pt=this->node_pt(j);
85 
86  // Boundary 1 is left (inflow) boundary
87  if (nod_pt->x(0)<0.226)
88  {
89  // If it's not on the upper or lower channel walls remove it
90  // from boundary 0
91  if ((nod_pt->x(1)<4.08)&&(nod_pt->x(1)>0.113))
92  {
93  this->remove_boundary_node(0,nod_pt);
94  }
95  this->add_boundary_node(1,nod_pt);
96  }
97 
98  // Boundary 2 is right (outflow) boundary
99  if (nod_pt->x(0)>8.28)
100  {
101  // If it's not on the upper or lower channel walls remove it
102  // from boundary 0
103  if ((nod_pt->x(1)<4.08)&&(nod_pt->x(1)>0.113))
104  {
105  this->remove_boundary_node(0,nod_pt);
106  }
107  this->add_boundary_node(2,nod_pt);
108  }
109  }
110  this->setup_boundary_element_info();
111  }
112 
113  /// Empty Destructor
114  virtual ~ElasticTriangleMesh() { }
115 
116 
117 };
118 
119 
120 
121 //==start_namespace==============================
122 /// Namespace for physical parameters
123 //==================================================
124 namespace Global_Physical_Variables
125 {
126 
127  /// Reynolds number
128  double Re=0.0;
129 
130  /// Pseudo-solid Poisson ratio
131  double Nu=0.3;
132 
133  /// Constitutive law used to determine the mesh deformation
134  ConstitutiveLaw *Constitutive_law_pt=0;
135 
136 } // end_of_namespace
137 
138 
139 
140 ////////////////////////////////////////////////////////////////////////
141 ////////////////////////////////////////////////////////////////////////
142 ////////////////////////////////////////////////////////////////////////
143 
144 
145 
146 //==start_of_problem_class============================================
147 /// Unstructured fluid problem
148 //====================================================================
149 template<class ELEMENT>
150 class UnstructuredFluidProblem : public Problem
151 {
152 
153 public:
154 
155  /// Constructor
157 
158  /// Destructor (empty)
160 
161  /// Access function for the fluid mesh
163  {
164  return Fluid_mesh_pt;
165  }
166 
167  /// Set the boundary conditions
168  void set_boundary_conditions();
169 
170  /// Doc the solution
171  void doc_solution(DocInfo& doc_info);
172 
173 private:
174 
175  /// Fluid mesh
177 
178 
179 }; // end_of_problem_class
180 
181 
182 //==start_constructor=====================================================
183 /// Constructor
184 //========================================================================
185 template<class ELEMENT>
187 {
188 
189  //Create fluid mesh
190  string fluid_node_file_name="fluid.fig.1.node";
191  string fluid_element_file_name="fluid.fig.1.ele";
192  string fluid_poly_file_name="fluid.fig.1.poly";
193  Fluid_mesh_pt = new ElasticTriangleMesh<ELEMENT>(fluid_node_file_name,
194  fluid_element_file_name,
195  fluid_poly_file_name);
196 
197  //Set the boundary conditions
198  this->set_boundary_conditions();
199 
200  // Add fluid mesh
201  add_sub_mesh(fluid_mesh_pt());
202 
203  // Build global mesh
204  build_global_mesh();
205 
206  // Setup equation numbering scheme
207  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
208 
209 } // end_of_constructor
210 
211 
212 
213 //============start_of_set_boundary_conditions=============================
214 /// Set the boundary conditions for the problem and document the boundary
215 /// nodes
216 //==========================================================================
217 template<class ELEMENT>
219 {
220  // Doc pinned nodes
221  std::ofstream solid_bc_file("pinned_solid_nodes.dat");
222  std::ofstream u_bc_file("pinned_u_nodes.dat");
223  std::ofstream v_bc_file("pinned_v_nodes.dat");
224 
225  // Set the boundary conditions for fluid problem: All nodes are
226  // free by default -- just pin the ones that have Dirichlet conditions
227  // here.
228  unsigned nbound=Fluid_mesh_pt->nboundary();
229  for(unsigned ibound=0;ibound<nbound;ibound++)
230  {
231  unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
232  for (unsigned inod=0;inod<num_nod;inod++)
233  {
234  // Get node
235  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
236 
237  // Pin velocity everywhere apart from outlet where we
238  // have parallel outflow
239  if (ibound!=2)
240  {
241  // Pin it
242  nod_pt->pin(0);
243  // Doc it as pinned
244  u_bc_file << nod_pt->x(0) << " " << nod_pt->x(1) << std::endl;
245  }
246  // Pin it
247  nod_pt->pin(1);
248  // Doc it as pinned
249  v_bc_file << nod_pt->x(0) << " " << nod_pt->x(1) << std::endl;
250 
251  // Pin pseudo-solid positions everywhere
252  for(unsigned i=0;i<2;i++)
253  {
254  // Pin the node
255  SolidNode* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
256  nod_pt->pin_position(i);
257 
258  // Doc it as pinned
259  solid_bc_file << nod_pt->x(i) << " ";
260  }
261  solid_bc_file << std::endl;
262  }
263  } // end loop over boundaries
264 
265  // Close
266  solid_bc_file.close();
267 
268  // Complete the build of all elements so they are fully functional
269  unsigned n_element = fluid_mesh_pt()->nelement();
270  for(unsigned e=0;e<n_element;e++)
271  {
272  // Upcast from GeneralisedElement to the present element
273  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(e));
274 
275  //Set the Reynolds number
276  el_pt->re_pt() = &Global_Physical_Variables::Re;
277 
278  // Set the constitutive law
279  el_pt->constitutive_law_pt() =
281  }
282 
283 
284  // Apply fluid boundary conditions: Poiseuille at inflow
285  //------------------------------------------------------
286 
287  // Find max. and min y-coordinate at inflow
288  unsigned ibound=1;
289  //Initialise y_min and y_max to y-coordinate of first node on boundary
290  double y_min=fluid_mesh_pt()->boundary_node_pt(ibound,0)->x(1);
291  double y_max=y_min;
292  //Loop over the rest of the nodes
293  unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
294  for (unsigned inod=1;inod<num_nod;inod++)
295  {
296  double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
297  if (y>y_max)
298  {
299  y_max=y;
300  }
301  if (y<y_min)
302  {
303  y_min=y;
304  }
305  }
306  double y_mid=0.5*(y_min+y_max);
307 
308  // Loop over all boundaries
309  for (unsigned ibound=0;ibound<fluid_mesh_pt()->nboundary();ibound++)
310  {
311  unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
312  for (unsigned inod=0;inod<num_nod;inod++)
313  {
314  // Parabolic inflow at the left boundary (boundary 1)
315  if (ibound==1)
316  {
317  double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
318  double veloc=1.5/(y_max-y_min)*
319  (y-y_min)*(y_max-y)/((y_mid-y_min)*(y_max-y_mid));
320  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,veloc);
321  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
322  }
323  // Zero flow elsewhere apart from x-velocity on outflow boundary
324  else
325  {
326  if(ibound!=2)
327  {
328  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,0.0);
329  }
330  fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
331  }
332  }
333  }
334 } // end_of_set_boundary_conditions
335 
336 
337 
338 
339 //==start_of_doc_solution=================================================
340 /// Doc the solution
341 //========================================================================
342 template<class ELEMENT>
344 {
345  ofstream some_file;
346  char filename[100];
347 
348  // Number of plot points
349  unsigned npts;
350  npts=5;
351 
352  // Output solution
353  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
354  doc_info.number());
355  some_file.open(filename);
356  fluid_mesh_pt()->output(some_file,npts);
357  some_file.close();
358 } //end_of_doc_solution
359 
360 
361 ////////////////////////////////////////////////////////////////////////
362 ////////////////////////////////////////////////////////////////////////
363 ////////////////////////////////////////////////////////////////////////
364 
365 
366 
367 
368 
369 
370 
371 //==start_of_main======================================================
372 /// Driver for unstructured fluid test problem
373 //=====================================================================
374 int main()
375 {
376  // Label for output
377  DocInfo doc_info;
378 
379  //Set the constitutive law for the pseudo-elasticity
381  new GeneralisedHookean(&Global_Physical_Variables::Nu);
382 
383  //Taylor--Hood formulation
384  {
385  // Set output directory
386  doc_info.set_directory("RESLT_TH");
387 
388  // Step number
389  doc_info.number()=0;
390 
391  // Build the problem with TTaylorHoodElements
392  UnstructuredFluidProblem<PseudoSolidNodeUpdateElement<
393  TTaylorHoodElement<2>,
394  TPVDElement<2,3> > > problem;
395 
396  // Output boundaries
397  problem.fluid_mesh_pt()->output_boundaries("RESLT_TH/boundaries.dat");
398 
399  // Output the initial guess for the solution
400  problem.doc_solution(doc_info);
401  doc_info.number()++;
402 
403  // Parameter study
404  double re_increment=5.0;
405  unsigned nstep=2; // 10;
406  for (unsigned i=0;i<nstep;i++)
407  {
408  // Solve the problem
409  problem.newton_solve();
410 
411  // Output the solution
412  problem.doc_solution(doc_info);
413  doc_info.number()++;
414 
415  // Bump up Re
416  Global_Physical_Variables::Re+=re_increment;
417  } //end_of_parameter_study
418  }
419 
420 
421 
422  //Crouzeix--Raviart formulation
423  {
424  // Set output directory
425  doc_info.set_directory("RESLT_CR");
426 
427  // Step number
428  doc_info.number()=0;
429 
430  // Reset Reynolds number
432 
433  // Build the problem with TCrouzeixRaviartElements
434  UnstructuredFluidProblem<PseudoSolidNodeUpdateElement<
435  TCrouzeixRaviartElement<2>,
436  TPVDBubbleEnrichedElement<2,3> > > problem;
437 
438  // Output boundaries
439  problem.fluid_mesh_pt()->output_boundaries("RESLT_CR/boundaries.dat");
440 
441  // Output the initial guess for the solution
442  problem.doc_solution(doc_info);
443  doc_info.number()++;
444 
445  // Parameter study
446  double re_increment=5.0;
447  unsigned nstep=2; // 10;
448  for (unsigned i=0;i<nstep;i++)
449  {
450  // Solve the problem
451  problem.newton_solve();
452 
453  // Output the solution
454  problem.doc_solution(doc_info);
455  doc_info.number()++;
456 
457  // Bump up Re
458  Global_Physical_Variables::Re+=re_increment;
459  }
460  }
461 
462 
463 } // end_of_main
464 
465 
466 
467 
468 
469 
470 
471 
472 
473 
Unstructured fluid problem.
void set_boundary_conditions()
Set the boundary conditions.
~UnstructuredFluidProblem()
Destructor (empty)
virtual ~ElasticTriangleMesh()
Empty Destructor.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
ElasticTriangleMesh< ELEMENT > *& fluid_mesh_pt()
Access function for the fluid mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
int main()
Driver for unstructured fluid test problem.
void identify_boundaries()
Function used to identify the domain boundaries.
ElasticTriangleMesh< ELEMENT > * Fluid_mesh_pt
Fluid mesh.
Triangle-based mesh upgraded to become a (pseudo-) solid mesh.
ElasticTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor:
double Nu
Pseudo-solid Poisson ratio.