driven_cavity.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-2023 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 //Driver for 2D rectangular driven cavity
27 
28 //Generic includes
29 #include "generic.h"
30 #include "navier_stokes.h"
31 
32 #include "meshes/simple_rectangular_quadmesh.h"
33 
34 
35 using namespace std;
36 
37 using namespace oomph;
38 
39 
40 //==start_of_namespace==============================
41 /// Namespace for physical parameters
42 //==================================================
44 {
45 
46  /// Reynolds number
47  double Re=100;
48 
49 } // end_of_namespace
50 
51 
52 
53 /// /////////////////////////////////////////////////////////////////////
54 /// /////////////////////////////////////////////////////////////////////
55 /// /////////////////////////////////////////////////////////////////////
56 
57 
58 
59 //==start_of_problem_class============================================
60 /// Driven cavity problem in rectangular domain
61 //====================================================================
62 template<class ELEMENT>
63 class RectangularDrivenCavityProblem : public Problem
64 {
65 
66 public:
67 
68 
69  /// Constructor
71 
72  /// Destructor (empty)
74 
75  /// Fix pressure in element e at pressure dof pdof and set to pvalue
76  void fix_pressure(const unsigned &e, const unsigned &pdof,
77  const double &pvalue)
78  {
79  //Cast to full element type and fix the pressure at that element
80  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
81  fix_pressure(pdof,pvalue);
82  } // end of fix_pressure
83 
84 
85  /// Update the after solve (empty)
87 
88 
89  /// Update the problem specs before solve.
90  /// Re-set velocity boundary conditions just to be on the safe side...
92  {
93  // Setup tangential flow along boundary 0:
94  unsigned ibound=0;
95  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
96  for (unsigned inod=0;inod<num_nod;inod++)
97  {
98  // Tangential flow
99  unsigned i=0;
100  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,1.0);
101  // No penetration
102  i=1;
103  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
104  }
105 
106  // Overwrite with no flow along the other boundaries
107  unsigned num_bound = mesh_pt()->nboundary();
108  for(unsigned ibound=1;ibound<num_bound;ibound++)
109  {
110  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
111  for (unsigned inod=0;inod<num_nod;inod++)
112  {
113  for (unsigned i=0;i<2;i++)
114  {
115  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(i,0.0);
116  }
117  }
118  }
119  } // end_of_actions_before_newton_solve
120 
121  // Access function for the specific mesh
122  SimpleRectangularQuadMesh<ELEMENT>* mesh_pt()
123  {
124  // Upcast from pointer to the Mesh base class to the specific
125  // element type that we're using here.
126  return dynamic_cast<SimpleRectangularQuadMesh<ELEMENT>*>(
127  Problem::mesh_pt());
128  }
129 
130 
131  /// Doc the solution
132  void doc_solution(DocInfo& doc_info);
133 
134 }; // end_of_problem_class
135 
136 
137 //==start_of_constructor==================================================
138 /// Constructor for RectangularDrivenCavity problem
139 //========================================================================
140 template<class ELEMENT>
142 {
143 
144  // Setup mesh
145 
146  // # of elements in x-direction
147  unsigned n_x=10;
148 
149  // # of elements in y-direction
150  unsigned n_y=10;
151 
152  // Domain length in x-direction
153  double l_x=1.0;
154 
155  // Domain length in y-direction
156  double l_y=1.0;
157 
158  // Build and assign mesh
159  Problem::mesh_pt() = new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
160 
161  // Set the boundary conditions for this problem: All nodes are
162  // free by default -- just pin the ones that have Dirichlet conditions
163  // here.
164  unsigned num_bound = mesh_pt()->nboundary();
165  for(unsigned ibound=0;ibound<num_bound;ibound++)
166  {
167  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
168  for (unsigned inod=0;inod<num_nod;inod++)
169  {
170  // Loop over values (u and v velocities)
171  for (unsigned i=0;i<2;i++)
172  {
173  mesh_pt()->boundary_node_pt(ibound,inod)->pin(i);
174  }
175  }
176  } // end loop over boundaries
177 
178  // Complete the build of all elements so they are fully functional
179 
180  //Find number of elements in mesh
181  unsigned n_element = mesh_pt()->nelement();
182 
183  // Loop over the elements to set up element-specific
184  // things that cannot be handled by constructor
185  for(unsigned e=0;e<n_element;e++)
186  {
187  // Upcast from GeneralisedElement to the present element
188  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
189 
190  //Set the Reynolds number
191  el_pt->re_pt() = &Global_Physical_Variables::Re;
192  } // end loop over elements
193 
194  // Now set the first pressure value in element 0 to 0.0
195  fix_pressure(0,0,0.0);
196 
197  // Setup equation numbering scheme
198  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
199 } // end_of_constructor
200 
201 
202 
203 //==start_of_doc_solution=================================================
204 /// Doc the solution
205 //========================================================================
206 template<class ELEMENT>
208 {
209  ofstream some_file;
210  char filename[100];
211 
212  // Number of plot points
213  unsigned npts;
214  npts=5;
215 
216  // Output solution
217  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
218  doc_info.number());
219  some_file.open(filename);
220  mesh_pt()->output(some_file,npts);
221  some_file.close();
222 } // end_of_doc_solution
223 
224 
225 
226 
227 
228 /// /////////////////////////////////////////////////////////////////////
229 /// /////////////////////////////////////////////////////////////////////
230 /// /////////////////////////////////////////////////////////////////////
231 
232 
233 
234 
235 
236 
237 
238 //==start_of_main======================================================
239 /// Driver for RectangularDrivenCavity test problem -- test drive
240 /// with two different types of element.
241 //=====================================================================
242 int main()
243 {
244 
245  // Set up doc info
246  // ---------------
247 
248  // Label for output
249  DocInfo doc_info;
250 
251  // Set output directory
252  doc_info.set_directory("RESLT");
253 
254  // Step number
255  doc_info.number()=0;
256 
257  // ---------------
258  // end of Set up doc info
259 
260  // Doing QCrouzeixRaviartElements
261  {
262 
263  // Build the problem with QCrouzeixRaviartElements
265  cout << "Doing QCrouzeixRaviartElement<2>" << std::endl;
266 
267  // Solve the problem
268  problem.newton_solve();
269 
270  // Outpt the solution
271  problem.doc_solution(doc_info);
272 
273  // Step number
274  doc_info.number()++;
275 
276  } // end of QCrouzeixRaviartElements
277 
278  // Doing QTaylorHoodElements
279  {
280 
281  // Build the problem with QTaylorHoodElements
283  cout << "Doing QTaylorHoodElement<2>" << std::endl;
284 
285  // Solve the problem
286  problem.newton_solve();
287 
288  // Outpt the solution
289  problem.doc_solution(doc_info);
290 
291  // Step number
292  doc_info.number()++;
293 
294  } // end of QTaylorHoodElements
295 
296 
297 } // end_of_main
298 
299 
300 
301 
302 
303 
304 
305 
306 
307 
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
~RectangularDrivenCavityProblem()
Destructor (empty)
SimpleRectangularQuadMesh< ELEMENT > * mesh_pt()
void actions_before_newton_solve()
Update the problem specs before solve. Re-set velocity boundary conditions just to be on the safe sid...
void doc_solution(DocInfo &doc_info)
Doc the solution.
void actions_after_newton_solve()
Update the after solve (empty)
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
RectangularDrivenCavityProblem()
Constructor.
int main()
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Namespace for physical parameters.
double Re
Reynolds number.