unstructured_three_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 code for a simple unstructured fluid problem using a mesh
31 // generated from an input file generated by the 3d mesh generator
32 // tetgen
33 
34 
35 //Generic routines
36 #include "generic.h"
37 #include "constitutive.h"
38 #include "navier_stokes.h"
39 
40 // Get the mesh
41 #include "meshes/tetgen_mesh.h"
42 
43 using namespace std;
44 using namespace oomph;
45 
46 
47 //////////////////////////////////////////////////////////////////
48 //////////////////////////////////////////////////////////////////
49 //////////////////////////////////////////////////////////////////
50 
51 
52 //=======start_namespace==========================================
53 /// Global variables
54 //================================================================
55 namespace Global_Parameters
56 {
57 
58  /// Default Reynolds number
59  double Re=100.0;
60 
61  /// Fluid pressure on inflow boundary
62  double P_in=0.5;
63 
64  /// Applied traction on fluid at the inflow boundary
65  void prescribed_inflow_traction(const double& t,
66  const Vector<double>& x,
67  const Vector<double>& n,
68  Vector<double>& traction)
69  {
70  traction[0]=0.0;
71  traction[1]=0.0;
72  traction[2]=P_in;
73  }
74 
75 
76  /// Fluid pressure on outflow boundary
77  double P_out=-0.5;
78 
79  /// Applied traction on fluid at the inflow boundary
80  void prescribed_outflow_traction(const double& t,
81  const Vector<double>& x,
82  const Vector<double>& n,
83  Vector<double>& traction)
84  {
85  traction[0]=0.0;
86  traction[1]=0.0;
87  traction[2]=-P_out;
88  }
89 
90 } //end namespace
91 
92 
93 
94 
95 
96 
97 //======start_problem_class===========================================
98 /// Unstructured fluid problem
99 //====================================================================
100 template<class ELEMENT>
101 class UnstructuredFluidProblem : public Problem
102 {
103 
104 public:
105 
106  /// Constructor:
108 
109  /// Destructor (empty)
111 
112  /// Doc the solution
113  void doc_solution(DocInfo& doc_info);
114 
115  /// Return total number of fluid inflow traction boundaries
117  {
118  return Inflow_boundary_id.size();
119  }
120 
121  /// Return total number of fluid outflow traction boundaries
123  {
124  return Outflow_boundary_id.size();
125  }
126 
127  /// Return total number of fluid outflow traction boundaries
129  {
130  return Inflow_boundary_id.size()+Outflow_boundary_id.size();
131  }
132 
133  //private:
134 
135  /// Create fluid traction elements at inflow
136  void create_fluid_traction_elements();
137 
138  /// Bulk fluid mesh
139  TetgenMesh<ELEMENT>* Fluid_mesh_pt;
140 
141  /// Meshes of fluid traction elements that apply pressure at in/outflow
142  Vector<Mesh*> Fluid_traction_mesh_pt;
143 
144  /// \short IDs of fluid mesh boundaries along which inflow boundary conditions
145  /// are applied
146  Vector<unsigned> Inflow_boundary_id;
147 
148  /// \short IDs of fluid mesh boundaries along which inflow boundary conditions
149  /// are applied
150  Vector<unsigned> Outflow_boundary_id;
151 
152 };
153 
154 
155 
156 //==========start_constructor=============================================
157 /// Constructor for unstructured 3D fluid problem
158 //========================================================================
159 template<class ELEMENT>
161 {
162 
163  //Create fluid bulk mesh, sub-dividing "corner" elements
164  string node_file_name="fsi_bifurcation_fluid.1.node";
165  string element_file_name="fsi_bifurcation_fluid.1.ele";
166  string face_file_name="fsi_bifurcation_fluid.1.face";
167  bool split_corner_elements=true;
168  Fluid_mesh_pt = new TetgenMesh<ELEMENT>(node_file_name,
169  element_file_name,
170  face_file_name,
171  split_corner_elements);
172 
173  // Find elements next to boundaries
174  //Fluid_mesh_pt->setup_boundary_element_info();
175 
176  // The following corresponds to the boundaries as specified by
177  // facets in the tetgen input:
178 
179  // Fluid mesh has one inflow boundary: Boundary 0
180  Inflow_boundary_id.resize(1);
181  Inflow_boundary_id[0]=0;
182 
183  // Fluid mesh has two outflow boundaries: Boundaries 1 and 2
184  Outflow_boundary_id.resize(2);
185  Outflow_boundary_id[0]=1;
186  Outflow_boundary_id[1]=2;
187 
188  // Apply BCs
189  //----------
190 
191  // Map to indicate which boundary has been done
192  std::map<unsigned,bool> done;
193 
194  // Loop over inflow/outflow boundaries to impose parallel flow
195  for (unsigned in_out=0;in_out<2;in_out++)
196  {
197  // Loop over in/outflow boundaries
198  unsigned n=nfluid_inflow_traction_boundary();
199  if (in_out==1) n=nfluid_outflow_traction_boundary();
200  for (unsigned i=0;i<n;i++)
201  {
202  // Get boundary ID
203  unsigned b=0;
204  if (in_out==0)
205  {
206  b=Inflow_boundary_id[i];
207  }
208  else
209  {
210  b=Outflow_boundary_id[i];
211  }
212 
213  // Number of nodes on that boundary
214  unsigned num_nod=Fluid_mesh_pt->nboundary_node(b);
215  for (unsigned inod=0;inod<num_nod;inod++)
216  {
217  // Get the node
218  Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(b,inod);
219 
220  // Pin transverse velocities
221  nod_pt->pin(0);
222  nod_pt->pin(1);
223  }
224 
225  // Done!
226  done[b]=true;
227  }
228 
229  } // done in and outflow
230 
231 
232 
233  // Loop over all fluid mesh boundaries and pin velocities
234  // of nodes that haven't been dealt with yet
235  unsigned nbound=Fluid_mesh_pt->nboundary();
236  for(unsigned b=0;b<nbound;b++)
237  {
238 
239  // Has the boundary been done yet?
240  if (!done[b])
241  {
242  unsigned num_nod=Fluid_mesh_pt->nboundary_node(b);
243  for (unsigned inod=0;inod<num_nod;inod++)
244  {
245  // Get node
246  Node* nod_pt= Fluid_mesh_pt->boundary_node_pt(b,inod);
247 
248  // Pin all velocities
249  nod_pt->pin(0);
250  nod_pt->pin(1);
251  nod_pt->pin(2);
252  }
253  }
254 
255  } // done no slip elsewhere
256 
257 
258  // Complete the build of the fluid elements so they are fully functional
259  //----------------------------------------------------------------------
260  unsigned n_element = Fluid_mesh_pt->nelement();
261  for(unsigned e=0;e<n_element;e++)
262  {
263 
264  // Upcast from GeneralisedElement to the present element
265  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
266 
267  //Set the Reynolds number
268  el_pt->re_pt() = &Global_Parameters::Re;
269 
270  }
271 
272 
273  // Create meshes of fluid traction elements at inflow/outflow
274  //-----------------------------------------------------------
275 
276  // Create the meshes
277  unsigned n=nfluid_traction_boundary();
278  Fluid_traction_mesh_pt.resize(n);
279  for (unsigned i=0;i<n;i++)
280  {
281  Fluid_traction_mesh_pt[i]=new Mesh;
282  }
283 
284  // Populate them with elements
285  create_fluid_traction_elements();
286 
287 
288  // Combine the lot
289  //----------------
290 
291  // Add sub meshes:
292 
293  // Fluid bulk mesh
294  add_sub_mesh(Fluid_mesh_pt);
295 
296  // The fluid traction meshes
297  n=nfluid_traction_boundary();
298  for (unsigned i=0;i<n;i++)
299  {
300  add_sub_mesh(Fluid_traction_mesh_pt[i]);
301  }
302 
303  // Build global mesh
304  build_global_mesh();
305 
306  // Setup equation numbering scheme
307  std::cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
308 
309 } // end constructor
310 
311 
312 
313 //============start_of_fluid_traction_elements==============================
314 /// Create fluid traction elements
315 //=======================================================================
316 template<class ELEMENT>
318 {
319 
320  // Counter for number of fluid traction meshes
321  unsigned count=0;
322 
323  // Loop over inflow/outflow boundaries
324  for (unsigned in_out=0;in_out<2;in_out++)
325  {
326  // Loop over boundaries with fluid traction elements
327  unsigned n=nfluid_inflow_traction_boundary();
328  if (in_out==1) n=nfluid_outflow_traction_boundary();
329  for (unsigned i=0;i<n;i++)
330  {
331  // Get boundary ID
332  unsigned b=0;
333  if (in_out==0)
334  {
335  b=Inflow_boundary_id[i];
336  }
337  else
338  {
339  b=Outflow_boundary_id[i];
340  }
341 
342  // How many bulk elements are adjacent to boundary b?
343  unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
344 
345  // Loop over the bulk elements adjacent to boundary b
346  for(unsigned e=0;e<n_element;e++)
347  {
348  // Get pointer to the bulk element that is adjacent to boundary b
349  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
350  Fluid_mesh_pt->boundary_element_pt(b,e));
351 
352  //What is the index of the face of the element e along boundary b
353  int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
354 
355  // Create new element
356  NavierStokesTractionElement<ELEMENT>* el_pt=
357  new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,
358  face_index);
359 
360  // Add it to the mesh
361  Fluid_traction_mesh_pt[count]->add_element_pt(el_pt);
362 
363  // Set the pointer to the prescribed traction function
364  if (in_out==0)
365  {
366  el_pt->traction_fct_pt() =
368  }
369  else
370  {
371  el_pt->traction_fct_pt() =
373  }
374  }
375  // Bump up counter
376  count++;
377  }
378  }
379 
380  } // end of create_traction_elements
381 
382 
383 
384 //========================================================================
385 /// Doc the solution
386 //========================================================================
387 template<class ELEMENT>
389 {
390 
391  ofstream some_file;
392  char filename[100];
393 
394  // Number of plot points
395  unsigned npts;
396  npts=5;
397 
398 
399  // Output fluid solution
400  sprintf(filename,"%s/fluid_soln%i.dat",doc_info.directory().c_str(),
401  doc_info.number());
402  some_file.open(filename);
403  Fluid_mesh_pt->output(some_file,npts);
404  some_file.close();
405 
406 }
407 
408 
409 
410 
411 
412 //=============start_main=================================================
413 /// Demonstrate how to solve an unstructured 3D fluids problem
414 //========================================================================
415 int main(int argc, char **argv)
416 {
417  // Store command line arguments
418  CommandLineArgs::setup(argc,argv);
419 
420  // Label for output
421  DocInfo doc_info;
422 
423  // Parameter study
424  double Re_increment=100.0;
425  unsigned nstep=4;
426  if (CommandLineArgs::Argc==2)
427  {
428  std::cout << "Validation -- only doing two steps" << std::endl;
429  nstep=2;
430  }
431 
432 
433  //Taylor--Hood
434  {
435  // Output directory
436  doc_info.set_directory("RESLT_TH");
437 
438  //Set up the problem
440 
441  //Output initial guess
442  problem.doc_solution(doc_info);
443  doc_info.number()++;
444 
445  // Parameter study: Crank up the pressure drop along the vessel
446  for (unsigned istep=0;istep<nstep;istep++)
447  {
448  // Solve the problem
449  problem.newton_solve();
450 
451  //Output solution
452  problem.doc_solution(doc_info);
453  doc_info.number()++;
454 
455  // Bump up Reynolds number (equivalent to increasing the imposed pressure
456  // drop)
457  Global_Parameters::Re+=Re_increment;
458  }
459  }
460 
461  //Crouzeix Raviart
462  {
463  //Reset to default Reynolds number
464  Global_Parameters::Re = 100.0;
465 
466  //Reset doc info number
467  doc_info.number()=0;
468 
469  // Output directory
470  doc_info.set_directory("RESLT_CR");
471 
472  //Set up the problem
474 
475  //Output initial guess
476  problem.doc_solution(doc_info);
477  doc_info.number()++;
478 
479  // Parameter study: Crank up the pressure drop along the vessel
480  for (unsigned istep=0;istep<nstep;istep++)
481  {
482  // Solve the problem
483  problem.newton_solve();
484 
485  //Output solution
486  problem.doc_solution(doc_info);
487  doc_info.number()++;
488 
489  // Bump up Reynolds number (equivalent to increasing the imposed pressure
490  // drop)
491  Global_Parameters::Re+=Re_increment;
492  }
493  }
494 
495 } // end_of_main
496 
497 
498 
499 
Unstructured fluid problem.
double P_in
Fluid pressure on inflow boundary.
void prescribed_inflow_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Applied traction on fluid at the inflow boundary.
TetgenMesh< ELEMENT > * Fluid_mesh_pt
Bulk fluid mesh.
Vector< unsigned > Outflow_boundary_id
IDs of fluid mesh boundaries along which inflow boundary conditions are applied.
~UnstructuredFluidProblem()
Destructor (empty)
unsigned nfluid_traction_boundary()
Return total number of fluid outflow traction boundaries.
void create_fluid_traction_elements()
Create fluid traction elements at inflow.
void prescribed_outflow_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Applied traction on fluid at the inflow boundary.
double P_out
Fluid pressure on outflow boundary.
unsigned nfluid_inflow_traction_boundary()
Return total number of fluid inflow traction boundaries.
int main(int argc, char **argv)
Demonstrate how to solve an unstructured 3D fluids problem.
void doc_solution(DocInfo &doc_info)
Doc the solution.
Vector< Mesh * > Fluid_traction_mesh_pt
Meshes of fluid traction elements that apply pressure at in/outflow.
unsigned nfluid_outflow_traction_boundary()
Return total number of fluid outflow traction boundaries.
double Re
Default Reynolds number.
Vector< unsigned > Inflow_boundary_id
IDs of fluid mesh boundaries along which inflow boundary conditions are applied.