spherical_cap_in_cylinder.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 refineable 2D Young Laplace problem on a circle sector
27 
28 //Generic routines
29 #include "generic.h"
30 
31 // The YoungLaplace equations
32 #include "young_laplace.h"
33 
34 
35 // The mesh
36 #include "meshes/quarter_circle_sector_mesh.h"
37 
38 
39 using namespace std;
40 using namespace oomph;
41 
42 
43 // Namespace (shared with other driver codes)
45 
46 
47 //====== start_of_problem_class=======================================
48 /// 2D RefineableYoungLaplace problem on a circle sector, discretised with
49 /// 2D QRefineableYoungLaplace elements. The specific type of element is
50 /// specified via the template parameter.
51 //====================================================================
52 template<class ELEMENT>
53 class RefineableYoungLaplaceProblem : public Problem
54 {
55 
56 public:
57 
58  /// Constructor:
60 
61  /// Destructor (empty)
63 
64  /// Update the problem specs before solve: Empty
66 
67  /// Update the problem after solve: Empty
69 
70  /// Increment problem parameters
72 
73  /// Doc the solution. DocInfo object stores flags/labels for where the
74  /// output gets written to and the trace file
75  void doc_solution(DocInfo& doc_info, ofstream& trace_file);
76 
77 private:
78 
79  /// Pointer to GeomObject that specifies the domain bondary
80  GeomObject* Boundary_pt;
81 
82  /// Pointer to the "bulk" mesh
83  RefineableQuarterCircleSectorMesh<ELEMENT>* Bulk_mesh_pt;
84 
85  /// Pointer to mesh containing the height control element
86  Mesh* Height_control_mesh_pt;
87 
88  /// Pointer to height control element
89  HeightControlElement* Height_control_element_pt;
90 
91  /// Node at which the height (displacement along spine) is controlled/doced
92  Node* Control_node_pt;
93 
94 }; // end of problem class
95 
96 
97 //=====start_of_constructor===============================================
98 /// Constructor for RefineableYoungLaplace problem
99 //========================================================================
100 template<class ELEMENT>
102 {
103 
104  // Setup dependent parameters in namespace
106 
107  // Setup bulk mesh
108  //----------------
109 
110  // Build geometric object that forms the curvilinear domain boundary:
111  // a unit circle
112 
113  // Create GeomObject
114  Boundary_pt=new Circle(0.0,0.0,1.0);
115 
116  // Start and end coordinates of curvilinear domain boundary on circle
117  double xi_lo=0.0;
118  double xi_hi=MathematicalConstants::Pi/2.0;
119 
120  // Now create the bulk mesh. Separating line between the two
121  // elements next to the curvilinear boundary is located half-way
122  // along the boundary.
123  double fract_mid=0.5;
124  Bulk_mesh_pt = new RefineableQuarterCircleSectorMesh<ELEMENT>(
125  Boundary_pt,xi_lo,fract_mid,xi_hi);
126 
127  // Create/set error estimator
128  Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
129 
130  // Set targets for spatial adaptivity
131  Bulk_mesh_pt->max_permitted_error()=1.0e-4;
132  Bulk_mesh_pt->min_permitted_error()=1.0e-6;
133 
134  // Add bulk mesh to the global mesh
135  add_sub_mesh(Bulk_mesh_pt);
136 
137 
138  // Prescribed height?
139  //-------------------
140 
141 
142  // Which element are we using for displacement control?
144 
145  // Choose the prescribed height element
146  ELEMENT* prescribed_height_element_pt= dynamic_cast<ELEMENT*>(
147  Bulk_mesh_pt->element_pt(GlobalParameters::Control_element));
148 
149  // ...and the associated control node (node 0 in that element)
150  // (we're storing this node even if there's no height-control, for
151  // output purposes...)
152  Control_node_pt= static_cast<Node*>(
153  prescribed_height_element_pt->node_pt(0));
154 
155  cout << "Controlling height at (x,y) : (" << Control_node_pt->x(0)
156  << "," << Control_node_pt->x(1) << ")" << endl;
157 
158  // If needed, create a height control element and store the
159  // pointer to the Kappa Data created by this object
160  Height_control_element_pt=0;
161  Height_control_mesh_pt=0;
163  {
164  Height_control_element_pt=new HeightControlElement(
165  Control_node_pt,&GlobalParameters::Controlled_height);
166 
167  GlobalParameters::Kappa_pt=Height_control_element_pt->kappa_pt();
168 
169  // Add to mesh
170  Height_control_mesh_pt = new Mesh;
171  Height_control_mesh_pt->add_element_pt(Height_control_element_pt);
172 
173  // Add height control mesh to the global mesh
174  add_sub_mesh(Height_control_mesh_pt);
175 
176  }
177  //...otherwise create a kappa data item from scratch and pin its
178  // single unknown value
179  else
180  {
181  GlobalParameters::Kappa_pt=new Data(1);
184  }
185 
186  // Build global mesh
187  //------------------
188  build_global_mesh();
189 
190 
191  // Boundary conditions
192  //--------------------
193 
194  // Set the boundary conditions for this problem: All nodes are
195  // free by default -- only need to pin the ones that have Dirichlet conditions
196  // here.
197  unsigned n_node = Bulk_mesh_pt->nboundary_node(1);
198  for (unsigned n=0;n<n_node;n++)
199  {
200  Bulk_mesh_pt->boundary_node_pt(1,n)->pin(0);
201  }
202 
203 
204  // Complete build of elements
205  //---------------------------
206 
207  // Complete the build of all elements so they are fully functional
208  unsigned n_bulk=Bulk_mesh_pt->nelement();
209  for(unsigned i=0;i<n_bulk;i++)
210  {
211  // Upcast from GeneralsedElement to the present element
212  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
213 
215  {
216  //Set the spine function pointers
217  el_pt->spine_base_fct_pt() = GlobalParameters::spine_base_function;
218  el_pt->spine_fct_pt() = GlobalParameters::spine_function;
219  }
220 
221  // Set the curvature data for the element
222  el_pt->set_kappa(GlobalParameters::Kappa_pt);
223 
224  }
225 
226  // Setup equation numbering scheme
227  cout <<"\nNumber of equations: " << assign_eqn_numbers() << endl;
228  cout << "\n********************************************\n" << endl;
229 
230 } // end of constructor
231 
232 
233 
234 //===============start_of_update_parameters==============================
235 /// Update (increase/decrease) parameters
236 //=======================================================================
237 template<class ELEMENT>
239 {
240 
243 
244  cout << "Solving for Prescribed Height Value = " ;
245  cout << GlobalParameters::Controlled_height << "\n" << endl;
246 
247 }
248 
249 
250 
251 //===============start_of_doc=============================================
252 /// Doc the solution: doc_info contains labels/output directory etc.
253 //========================================================================
254 template<class ELEMENT>
256  ofstream& trace_file)
257 {
258 
259  // Output kappa vs height
260  //-----------------------
261  trace_file << -1.0*GlobalParameters::Kappa_pt->value(0) << " ";
262  trace_file << GlobalParameters::get_exact_kappa() << " ";
263  trace_file << Control_node_pt->value(0) ;
264  trace_file << endl;
265 
266 
267  // Number of plot points: npts x npts
268  unsigned npts=5;
269 
270  // Output full solution
271  //---------------------
272  ofstream some_file;
273  char filename[100];
274  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
275  doc_info.number());
276  some_file.open(filename);
277  Bulk_mesh_pt->output(some_file,npts);
278  some_file.close();
279 
280 }
281 
282 
283 
284 //===== start_of_main=====================================================
285 /// Driver code for 2D RefineableYoungLaplace problem. Input arguments: none
286 /// (for validation) or number of steps.
287 //========================================================================
288 int main(int argc, char* argv[])
289 {
290 
291  // Store command line arguments
292  CommandLineArgs::setup(argc,argv);
293 
294  // No command line args: Running with limited number of steps
295  if (CommandLineArgs::Argc==1)
296  {
297  std::cout
298  << "Running with limited number of steps for validation"
299  << std::endl;
300 
301  // Number of steps
303  }
304  else
305  {
306  // Number of steps
307  GlobalParameters::Nsteps=atoi(argv[1]);
308  }
309  // Create label for output
310  //------------------------
311  DocInfo doc_info;
312 
313  // Set outputs
314  //------------
315 
316  // Trace file
317  ofstream trace_file;
318 
319 // Set output directory
320  doc_info.set_directory("RESLT_adapt_pinned_spherical_cap_in_cylinder");
321 
322 //Open a trace file
323  char filename[100];
324  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
325  trace_file.open(filename);
326 
327  trace_file
328  << "VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{exact}\",\"h\""
329  << std::endl;
330  trace_file << "ZONE" << std::endl;
331 
332 
333  // Set case
335 
336  // Run with spines
338 
339 
340  //Set up the problem
341  //------------------
342 
343  // Create the problem with 2D nine-node elements from the
344  // RefineableQYoungLaplaceElement family.
346 
347  //Output the solution
348  problem.doc_solution(doc_info,trace_file);
349 
350  //Increment counter for solutions
351  doc_info.number()++;
352 
353  // Parameter incrementation
354  //-------------------------
355 
356  // Loop over steps
357  for (unsigned istep=0;istep<GlobalParameters::Nsteps;istep++)
358  {
359  // Solve the problem
360  unsigned max_adapt=1;
361  problem.newton_solve(max_adapt);
362 
363  //Output the solution
364  problem.doc_solution(doc_info,trace_file);
365 
366  //Increment counter for solutions
367  doc_info.number()++;
368 
369  // Increase the parameters
370  problem.increment_parameters();
371  }
372 
373  // Close output file
374  trace_file.close();
375 
376 
377 
378 
379 
380 
381 
382 } //end of main
383 
384 
2D RefineableYoungLaplace problem on rectangular domain, discretised with 2D QRefineableYoungLaplace ...
~RefineableYoungLaplaceProblem()
Destructor (empty)
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...
void actions_after_solve()
Update the problem after solve: Empty.
GeomObject * Boundary_pt
Pointer to GeomObject that specifies the domain bondary.
RefineableYoungLaplaceProblem()
Constructor:
RefineableQuarterCircleSectorMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void increment_parameters()
Increment problem parameters.
void actions_before_solve()
Update the problem specs before solve: Empty.
double Controlled_height
Height control value.
Definition: barrel.cc:51
unsigned Control_element
Number of element in bulk mesh at which height control is applied. Initialise to 0 – will be overwrit...
double get_exact_kappa()
Exact kappa.
Definition: barrel.cc:54
bool Use_height_control
Use height control (true) or not (false)?
double Controlled_height_increment
Increment for height control.
bool Use_spines
Use spines (true) or not (false)
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
unsigned Nsteps
Number of steps.
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 Kappa_initial
Initial value for kappa.
int Case
What case are we considering: Choose one from the enumeration Cases.
void setup_dependent_parameters_and_sanity_check()
Setup dependent parameters and perform sanity check.
int main(int argc, char *argv[])
Driver code for 2D RefineableYoungLaplace problem. Input arguments: none (for validation) or number o...