barrel.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 code for a 2D YoungLaplace problem
27 
28 // Generic routines
29 #include "generic.h"
30 
31 // The YoungLaplace equations
32 #include "young_laplace.h"
33 
34 // The mesh
35 #include "meshes/simple_rectangular_quadmesh.h"
36 
37 
38 using namespace std;
39 using namespace oomph;
40 
41 //===== start_of_namespace========================================
42 /// Namespace for "global" problem parameters
43 //================================================================
45 {
46 
47  // Displacement control:
48  //----------------------
49 
50  /// Height control value
51  double Controlled_height = 0.0;
52 
53  /// Exact kappa
54  double get_exact_kappa()
55  {
56 
57  // Mean curvature of barrel-shaped meniscus
58  return 2.0*Controlled_height/
60 
61  } //end exact kappa
62 
63 
64  // Spine basis
65  //------------
66 
67  /// Spine basis: The position vector to the basis of the spine
68  /// as a function of the two coordinates x_1 and x_2, and its
69  /// derivatives w.r.t. to these coordinates.
70  /// dspine_B[i][j] = d spine_B[j] / dx_i
71  /// Spines start in the (x_1,x_2) plane at (x_1,x_2).
72  void spine_base_function(const Vector<double>& x,
73  Vector<double>& spine_B,
74  Vector< Vector<double> >& dspine_B)
75  {
76 
77  // Bspines and derivatives
78  spine_B[0] = x[0];
79  spine_B[1] = x[1];
80  spine_B[2] = 0.0 ;
81  dspine_B[0][0] = 1.0 ;
82  dspine_B[1][0] = 0.0 ;
83  dspine_B[0][1] = 0.0 ;
84  dspine_B[1][1] = 1.0 ;
85  dspine_B[0][2] = 0.0 ;
86  dspine_B[1][2] = 0.0 ;
87 
88  } // End of bspine functions
89 
90 
91 
92  // Spines rotate in the y-direction
93  //---------------------------------
94 
95  /// Min. spine angle against horizontal plane
96  double Alpha_min = MathematicalConstants::Pi/2.0*1.5;
97 
98  /// Max. spine angle against horizontal plane
99  double Alpha_max = MathematicalConstants::Pi/2.0*0.5;
100 
101  /// Spine: The spine vector field as a function of the two
102  /// coordinates x_1 and x_2, and its derivatives w.r.t. to these coordinates:
103  /// dspine[i][j] = d spine[j] / dx_i
104  void spine_function(const Vector<double>& x,
105  Vector<double>& spine,
106  Vector< Vector<double> >& dspine)
107  {
108 
109  /// Spines (and derivatives) are independent of x[0] and rotate
110  /// in the x[1]-direction
111  spine[0]=0.0;
112  dspine[0][0]=0.0;
113  dspine[1][0]=0.0;
114 
115  spine[1]=cos(Alpha_min+(Alpha_max-Alpha_min)*x[1]);
116  dspine[0][1]=0.0;
117  dspine[1][1]=-sin(Alpha_min+(Alpha_max-Alpha_min)*x[1])
118  *(Alpha_max-Alpha_min);
119 
120  spine[2]=sin(Alpha_min+(Alpha_max-Alpha_min)*x[1]);
121  dspine[0][2]=0.0;
122  dspine[1][2]=cos(Alpha_min+(Alpha_max-Alpha_min)*x[1])
123  *(Alpha_max-Alpha_min);
124 
125  } // End spine function
126 
127 
128 } // end of namespace
129 
130 
131 
132 
133 //====== start_of_problem_class=======================================
134 /// 2D YoungLaplace problem on rectangular domain, discretised with
135 /// 2D QYoungLaplace elements. The specific type of element is
136 /// specified via the template parameter.
137 //====================================================================
138 template<class ELEMENT>
139 class YoungLaplaceProblem : public Problem
140 {
141 
142 public:
143 
144  /// Constructor:
146 
147  /// Destructor (empty)
149 
150  /// Update the problem before solve
152  {
153  // This only has an effect if displacement control is disabled
154  double new_kappa=Kappa_pt->value(0)-0.5;
155  Kappa_pt->set_value(0,new_kappa);
156  }
157 
158  /// Update the problem after solve: Empty
160 
161  /// Doc the solution. DocInfo object stores flags/labels for where the
162  /// output gets written to and the trace file
163  void doc_solution(DocInfo& doc_info, ofstream& trace_file);
164 
165 private:
166 
167  /// Node at which the height (displacement along spine) is controlled/doced
169 
170  /// Pointer to Data object that stores the prescribed curvature
171  Data* Kappa_pt;
172 
173 }; // end of problem class
174 
175 
176 //=====start_of_constructor===============================================
177 /// Constructor for YoungLaplace problem
178 //========================================================================
179 template<class ELEMENT>
181 {
182 
183  // Setup mesh
184  //-----------
185 
186  // # of elements in x-direction
187  unsigned n_x=8;
188 
189  // # of elements in y-direction
190  unsigned n_y=8;
191 
192  // Domain length in x-direction
193  double l_x=1.0;
194 
195  // Domain length in y-direction
196  double l_y=1.0;
197 
198  // Build and assign mesh
199  Problem::mesh_pt()=new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
200 
201 
202  // Check that we've got an even number of elements otherwise
203  // out counting doesn't work...
204  if ((n_x%2!=0)||(n_y%2!=0))
205  {
206  cout << "n_x n_y should be even" << endl;
207  abort();
208  }
209 
210  // This is the element that contains the central node:
211  ELEMENT* prescribed_height_element_pt= dynamic_cast<ELEMENT*>(
212  mesh_pt()->element_pt(n_y*n_x/2+n_x/2));
213 
214  // The central node is node 0 in that element
215  Control_node_pt= static_cast<Node*>(prescribed_height_element_pt->node_pt(0));
216 
217  std::cout << "Controlling height at (x,y) : (" << Control_node_pt->x(0)
218  << "," << Control_node_pt->x(1) << ")" << "\n" << endl;
219 
220 
221  // Create a height control element
222  HeightControlElement* height_control_element_pt=new HeightControlElement(
223  Control_node_pt,&GlobalParameters::Controlled_height);
224 
225  // Store pointer to kappa data
226  Kappa_pt=height_control_element_pt->kappa_pt();
227 
228 
229  // Comment out the previous two commands and uncomment the following two
230  // to prescribe the pressure drop (the curvature) directly
231  //Kappa_pt=new Data(1);
232  //Kappa_pt->pin(0);
233 
234 
235  // Boundary conditions
236  //--------------------
237 
238  // Set the boundary conditions for this problem: All nodes are
239  // free by default -- only need to pin the ones that have Dirichlet conditions
240  // here.
241  unsigned n_bound = mesh_pt()->nboundary();
242  for(unsigned b=0;b<n_bound;b++)
243  {
244 
245  // Pin meniscus displacement at all nodes boundaries 0 and 2
246  if ((b==0)||(b==2))
247  {
248  unsigned n_node = mesh_pt()->nboundary_node(b);
249  for (unsigned n=0;n<n_node;n++)
250  {
251  mesh_pt()->boundary_node_pt(b,n)->pin(0);
252  }
253  }
254 
255  } // end bc
256 
257  // Complete build of elements
258  //---------------------------
259 
260  // Complete the build of all elements so they are fully functional
261  unsigned nelement = mesh_pt()->nelement();
262  for(unsigned i=0;i<nelement;i++)
263  {
264  // Upcast from GeneralsedElement to YoungLaplace element
265  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
266 
267  //Set the spine function pointers
268  el_pt->spine_base_fct_pt() = GlobalParameters::spine_base_function;
269  el_pt->spine_fct_pt() = GlobalParameters::spine_function;
270 
271  // Set the curvature data for the element
272  el_pt->set_kappa(Kappa_pt);
273  }
274 
275  // Add the height control element to mesh (comment this out
276  // if you're not using displacement control)
277  mesh_pt()->add_element_pt(height_control_element_pt);
278 
279  // Setup equation numbering scheme
280  cout <<"\nNumber of equations: " << assign_eqn_numbers() << endl;
281 
282 } // end of constructor
283 
284 
285 
286 
287 //===============start_of_doc=============================================
288 /// Doc the solution: doc_info contains labels/output directory etc.
289 //========================================================================
290 template<class ELEMENT>
292  ofstream& trace_file)
293 {
294 
295  // Output kappa vs height of the apex
296  //------------------------------------
297  trace_file << -1.0*Kappa_pt->value(0) << " ";
298  trace_file << GlobalParameters::get_exact_kappa() << " ";
299  trace_file << Control_node_pt->value(0) ;
300  trace_file << endl;
301 
302  // Number of plot points: npts x npts
303  unsigned npts=5;
304 
305  // Output full solution
306  //---------------------
307  ofstream some_file;
308  char filename[100];
309  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
310  doc_info.number());
311  some_file.open(filename);
312  mesh_pt()->output(some_file,npts);
313  some_file.close();
314 
315 } // end of doc
316 
317 
318 //===================start_of_main========================================
319 /// Driver code
320 //========================================================================
321 int main()
322 {
323 
324  // Create label for output
325  DocInfo doc_info;
326 
327  // Set output directory
328  doc_info.set_directory("RESLT");
329 
330  //Open a trace file
331  ofstream trace_file;
332  char filename[100];
333  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
334  trace_file.open(filename);
335 
336  // Write kappa, exact kappa and height values
337  trace_file
338  << "VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{ex}\",\"h\""
339  << std::endl;
340  trace_file << "ZONE" << std::endl;
341 
342  // Create the problem
343  //-------------------
344 
345  // Create the problem with 2D nine-node elements from the
346  // QYoungLaplaceElement family.
348 
349  //Output the solution
350  problem.doc_solution(doc_info,trace_file);
351 
352  //Increment counter for solutions
353  doc_info.number()++;
354 
355 
356  // Parameter incrementation
357  //-------------------------
358  double increment=0.1;
359 
360  // Loop over steps
361  unsigned nstep=2; // 10;
362  for (unsigned istep=0;istep<nstep;istep++)
363  {
364 
365  // Increment prescribed height value
367 
368  // Solve the problem
369  problem.newton_solve();
370 
371  //Output the solution
372  problem.doc_solution(doc_info,trace_file);
373 
374  //Increment counter for solutions
375  doc_info.number()++;
376 
377  }
378 
379  // Close output file
380  trace_file.close();
381 
382 } // end of main
383 
384 
385 
386 
387 
388 
int main()
Driver code.
Definition: barrel.cc:321
2D YoungLaplace problem on rectangular domain, discretised with 2D QYoungLaplace elements....
Definition: barrel.cc:140
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to and the tra...
Definition: barrel.cc:291
Node * Control_node_pt
Node at which the height (displacement along spine) is controlled/doced.
Definition: barrel.cc:168
YoungLaplaceProblem()
Constructor:
Definition: barrel.cc:180
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
Definition: barrel.cc:171
void actions_after_newton_solve()
Update the problem after solve: Empty.
Definition: barrel.cc:159
void actions_before_newton_solve()
Update the problem before solve.
Definition: barrel.cc:151
~YoungLaplaceProblem()
Destructor (empty)
Definition: barrel.cc:148
Namespace for "global" problem parameters.
Definition: barrel.cc:45
double Alpha_max
Max. spine angle against horizontal plane.
Definition: barrel.cc:99
double Controlled_height
Height control value.
Definition: barrel.cc:51
double get_exact_kappa()
Exact kappa.
Definition: barrel.cc:54
void spine_function(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double > > &dspine)
Spine: The spine vector field as a function of the two coordinates x_1 and x_2, and its derivatives w...
Definition: barrel.cc:104
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
void spine_base_function(const Vector< double > &x, Vector< double > &spine_B, Vector< Vector< double > > &dspine_B)
Spine basis: The position vector to the basis of the spine as a function of the two coordinates x_1 a...
Definition: barrel.cc:72
double Alpha_min
Min. spine angle against horizontal plane.
Definition: barrel.cc:96