two_d_unsteady_heat_adapt.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 adaptive unsteady heat problem
27 
28 //Generic routines
29 #include "generic.h"
30 
31 // The unsteady heat equations
32 #include "unsteady_heat.h"
33 
34 // Mesh
35 #include "meshes/quarter_circle_sector_mesh.h"
36 
37 using namespace std;
38 
39 using namespace oomph;
40 
41 using namespace MathematicalConstants;
42 
43 /// ////////////////////////////////////////////////////////////////////
44 /// ////////////////////////////////////////////////////////////////////
45 /// /////////////////////////////////////////////////////////////////////
46 
47 
48 //============start_of_MyUnitCircle====================================
49 /// Unit circle as GeomObject
50 /// \f[ x = \cos(\xi) \f]
51 /// \f[ y = \sin(\xi) \f]
52 //=====================================================================
53 class MyUnitCircle : public GeomObject
54 {
55 
56 public:
57 
58  /// Constructor: The circle is a 1D object
59  /// (i.e. it's parametrised by one intrinsic coordinate) in 2D space.
60  /// Pass these arguments to the constructor of the GeomObject base class.
61  MyUnitCircle() : GeomObject(1,2) {}
62 
63  /// Destructor: Empty
64  virtual ~MyUnitCircle() {}
65 
66  /// Current position vector to material point at
67  /// Lagrangian coordinate xi
68  void position(const Vector<double>& xi, Vector<double>& r) const
69  {
70  // Position vector
71  r[0] = cos(xi[0]);
72  r[1] = sin(xi[0]);
73  } // end of position(...)
74 
75 
76  /// Parametrised position on object: r(xi). Evaluated at
77  /// previous time level. t=0: current time; t>0: previous
78  /// time level. Circle is fixed -- simply call the steady version.
79  void position(const unsigned& t, const Vector<double>& xi,
80  Vector<double>& r) const
81  {
82  // Call steady version
83  position(xi,r);
84  } // end of position(...)
85 
86 
87 }; // end of MyUnitCircle
88 
89 
90 /// ////////////////////////////////////////////////////////////////////
91 /// ////////////////////////////////////////////////////////////////////
92 /// /////////////////////////////////////////////////////////////////////
93 
94 
95 //======start_of_TanhSolnForUnsteadyHeat==============================
96 /// Namespace for exact solution of unsteady heat equation
97 /// with sharp step
98 //====================================================================
100 {
101 
102  /// Parameter for steepness of step
103  double Alpha;
104 
105  /// Parameter for amplitude of step translation
106  double Beta;
107 
108  /// Parameter for timescale of step translation
109  double Gamma;
110 
111  /// Parameter for angle of step
112  double TanPhi;
113 
114  /// Position of step (x-axis intercept)
115  double step_position(const double& time)
116  {
117  return Beta*tanh(Gamma*cos(0.2E1*MathematicalConstants::Pi*time));
118  }
119 
120  /// Exact solution as a Vector
121  void get_exact_u(const double& time, const Vector<double>& x,
122  Vector<double>& u)
123  {
124  double X=x[0];
125  double Y=x[1];
126  u[0]=tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
127  MathematicalConstants::Pi*time)))-Y));
128  }
129 
130  /// Exact solution as a scalar
131  void get_exact_u(const double& time, const Vector<double>& x, double& u)
132  {
133  double X=x[0];
134  double Y=x[1];
135  u=tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
136  MathematicalConstants::Pi*time)))-Y));
137  }
138 
139  /// Source function to make it an exact solution
140  void get_source(const double& time, const Vector<double>& x, double& source)
141  {
142  double X=x[0];
143  double Y=x[1];
144  source = -2.0*tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
145 MathematicalConstants::Pi*time)))-Y))*(1.0-pow(tanh(0.1E1-Alpha*(TanPhi*(X-
146 Beta*tanh(Gamma*cos(0.2E1*MathematicalConstants::Pi*time)))-Y)),2.0))*Alpha*
147 Alpha*TanPhi*TanPhi-2.0*tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
148 MathematicalConstants::Pi*time)))-Y))*(1.0-pow(tanh(0.1E1-Alpha*(TanPhi*(X-
149 Beta*tanh(Gamma*cos(0.2E1*MathematicalConstants::Pi*time)))-Y)),2.0))*Alpha*
150 Alpha+0.2E1*(1.0-pow(tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
151 MathematicalConstants::Pi*time)))-Y)),2.0))*Alpha*TanPhi*Beta*(1.0-pow(tanh(
152 Gamma*cos(0.2E1*MathematicalConstants::Pi*time)),2.0))*Gamma*sin(0.2E1*
153 MathematicalConstants::Pi*time)*MathematicalConstants::Pi;
154  }
155 
156  /// Flux required by the exact solution on a boundary on which y is fixed
157  void prescribed_flux_on_fixed_y_boundary(const double& time,
158  const Vector<double>& x,
159  double& flux)
160  {
161  double X=x[0];
162  double Y=x[1];
163 
164  //The outer unit normal to the boundary is (0,-1)
165  double NX = 0.0;
166  double NY = -1.0;
167 
168  //The flux in terms of the normal is
169  flux = -(1.0-pow(tanh(0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*
170 MathematicalConstants::Pi*time)))-Y)),2.0))*Alpha*TanPhi*NX+(1.0-pow(tanh(
171 0.1E1-Alpha*(TanPhi*(X-Beta*tanh(Gamma*cos(0.2E1*MathematicalConstants::Pi*
172 time)))-Y)),2.0))*Alpha*NY;
173  }
174 
175 } // end of TanhSolnForUnsteadyHeat
176 
177 
178 /// /////////////////////////////////////////////////////////////////////
179 /// /////////////////////////////////////////////////////////////////////
180 /// /////////////////////////////////////////////////////////////////////
181 
182 
183 //=====start_of_problem_class=========================================
184 /// Unsteady heat problem in quarter circle domain.
185 //====================================================================
186 template<class ELEMENT>
187 class RefineableUnsteadyHeatProblem : public Problem
188 {
189 
190 public:
191 
192  /// Constructor: Pass pointer to source function
193  RefineableUnsteadyHeatProblem(UnsteadyHeatEquations<2>::
194  UnsteadyHeatSourceFctPt source_fct_pt);
195 
196  /// Destructor: Close trace file
198 
199  /// Update the problem specs after solve (empty)
201 
202  /// Update the problem specs before solve (empty)
204 
205  /// Update the problem specs after timestep (empty)
207 
208  /// Update the problem specs before next timestep:
209  /// Set Dirchlet boundary conditions from exact solution.
210  void actions_before_implicit_timestep();
211 
212  /// Actions before adapt: Wipe the mesh of prescribed flux elements
213  void actions_before_adapt();
214 
215  /// Actions after adapt: Rebuild the mesh of prescribed flux elements
216  void actions_after_adapt();
217 
218  /// Set initial condition (incl previous timesteps) according
219  /// to specified function. Note that his overloads the virtual
220  /// function in the Problem base class and is therefore executed
221  /// automatically to re-assign the initial conditions during the
222  /// spatially adaptive solution at the first timestep.
223  void set_initial_condition();
224 
225  /// Create UnsteadyHeat flux elements on boundary b of the Mesh pointed
226  /// to by bulk_mesh_pt and add them to the Mesh object pointed to by
227  /// surface_mesh_pt
228  void create_flux_elements(const unsigned &b, Mesh* const &bulk_mesh_pt,
229  Mesh* const &surface_mesh_pt);
230 
231  /// Delete UnsteadyHeat flux elements and wipe the surface mesh
232  void delete_flux_elements(Mesh* const &surface_mesh_pt);
233 
234  /// Doc the solution
235  void doc_solution();
236 
237  /// Dump problem data to allow for later restart
238  void dump_it(ofstream& dump_file);
239 
240  /// Read problem data for restart
241  void restart(ifstream& restart_file);
242 
243  /// Pointer to bulk mesh
244  RefineableQuarterCircleSectorMesh<ELEMENT>* bulk_mesh_pt()
245  {
246  return Bulk_mesh_pt;
247  }
248 
249 
250 private:
251 
252  /// Pointer to GeomObject that specifies the domain bondary
253  GeomObject* Boundary_pt;
254 
255  /// Pointer to source function
256  UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt Source_fct_pt;
257 
258  /// Pointer to the "bulk" mesh
259  RefineableQuarterCircleSectorMesh<ELEMENT>* Bulk_mesh_pt;
260 
261  /// Pointer to the "surface" mesh
263 
264  /// Pointer to central node (exists at all refinement levels) for doc
265  Node* Doc_node_pt;
266 
267  /// Doc info object
268  DocInfo Doc_info;
269 
270  /// Trace file
271  ofstream Trace_file;
272 
273 }; // end of problem_class
274 
275 //========start_of_constructor============================================
276 /// Constructor for UnsteadyHeat problem in quarter circle domain.
277 /// Pass pointer to source function.
278 //========================================================================
279 template<class ELEMENT>
281  UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt) :
282  Source_fct_pt(source_fct_pt)
283 {
284 
285 
286  // Setup labels for output
287  //------------------------
288 
289  // Output directory
290  Doc_info.set_directory("RESLT");
291 
292  // Output number
293  Doc_info.number()=0;
294 
295  // Open trace file
296  char filename[100];
297  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
298  Trace_file.open(filename);
299 
300  Trace_file << "VARIABLES=\"time t\",\"u<SUB>FE</SUB>\",\"u<SUB>exact</SUB>\","
301  << "\"A\","
302  << "\"X<SUB>step</SUB>\","
303  << "\"N<SUB>element</SUB>\","
304  << "\"N<SUB>refined</SUB>\","
305  << "\"N<SUB>unrefined</SUB>\","
306  << "\"norm of error\","
307  << "\"norm of solution\""
308  << std::endl;
309 
310 
311  // Setup parameters for tanh solution
312  // ----------------------------------
313 
314  // Steepness of step
316 
317  // Orientation of step
319 
320  // Amplitude for movement of step
322 
323  // Parameter for time-dependence of step movement
325 
326 
327  //Allocate the timestepper -- This constructs the time object as well
328  add_time_stepper_pt(new BDF<2>());
329 
330  // Setup mesh
331  //-----------
332 
333  // Build geometric object that forms the curvilinear domain boundary:
334  // a unit circle
335 
336  // Create GeomObject
338 
339  // Start and end coordinates of curvilinear domain boundary on circle
340  double xi_lo=0.0;
341  double xi_hi=MathematicalConstants::Pi/2.0;
342 
343  // Now create the bulk mesh. Separating line between the two
344  // elements next to the curvilinear boundary is located half-way
345  // along the boundary.
346  double fract_mid=0.5;
347  Bulk_mesh_pt = new RefineableQuarterCircleSectorMesh<ELEMENT>(
348  Boundary_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());
349 
350  // Create the surface mesh as an empty mesh
351  Surface_mesh_pt=new Mesh;
352 
353  // Create prescribed-flux elements from all elements that are
354  // adjacent to boundary 0 (the horizontal lower boundary), and add them
355  // to the (so far empty) surface mesh.
357 
358  // Add the two sub meshes to the problem
359  add_sub_mesh(Bulk_mesh_pt);
360  add_sub_mesh(Surface_mesh_pt);
361 
362  // Combine all submeshes into a single global Mesh
363  build_global_mesh();
364 
365  // Set error estimator for bulk mesh
366  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
367  Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
368 
369  // Set the boundary conditions for this problem: All nodes are
370  // free by default -- just pin the ones that have Dirichlet conditions
371  // here.
372  unsigned n_bound = Bulk_mesh_pt->nboundary();
373  for(unsigned b=0;b<n_bound;b++)
374  {
375  // Leave nodes on boundary 0 free -- this is where we apply the flux
376  // boundary condition
377  if (b!=0)
378  {
379  unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
380  for (unsigned n=0;n<n_node;n++)
381  {
382  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
383  }
384  }
385  }
386 
387  // Extract pointer to the central node (this exists at all refinement levels)
388  // for doc of solution
389  FiniteElement* el0_pt=Bulk_mesh_pt->finite_element_pt(0);
390  unsigned nnod=el0_pt->nnode();
391  Doc_node_pt=el0_pt->node_pt(nnod-1);
392 
393 
394  // Complete the build of all elements so they are fully functional
395  //----------------------------------------------------------------
396 
397  // Find number of elements in mesh
398  unsigned n_element = Bulk_mesh_pt->nelement();
399 
400  // Loop over the elements to set up element-specific
401  // things that cannot be handled by constructor
402  for(unsigned i=0;i<n_element;i++)
403  {
404  // Upcast from FiniteElement to the present element
405  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
406 
407  //Set the source function pointer
408  el_pt->source_fct_pt() = Source_fct_pt;
409  }
410 
411  // Loop over the flux elements to pass pointer to prescribed flux function
412  n_element=Surface_mesh_pt->nelement();
413  for(unsigned e=0;e<n_element;e++)
414  {
415  // Upcast from GeneralisedElement to UnsteadyHeat flux element
416  UnsteadyHeatFluxElement<ELEMENT> *el_pt =
417  dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*>(
418  Surface_mesh_pt->element_pt(e));
419 
420  // Set the pointer to the prescribed flux function
421  el_pt->flux_fct_pt() =
423  }
424 
425  // Do equation numbering
426  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
427 
428 } // end of constructor
429 
430 
431 //======start_of_destructor===============================================
432 /// Destructor: Close trace file
433 //========================================================================
434 template<class ELEMENT>
436 {
437  // Close trace file
438  Trace_file.close();
439 } // end of destructor
440 
441 
442 //=========start of actions_before_implicit_timestep===============================
443 /// Actions before timestep: Set the boundary conditions
444 /// for the current time.
445 //========================================================================
446 template<class ELEMENT>
448 {
449  // Get current time
450  double time=time_pt()->time();
451 
452  //Loop over all boundaries
453  unsigned num_bound = Bulk_mesh_pt->nboundary();
454  for(unsigned b=0;b<num_bound;b++)
455  {
456  // Exclude flux boundary
457  if (b!=0)
458  {
459  // Loop over the nodes on boundary
460  unsigned num_nod=Bulk_mesh_pt->nboundary_node(b);
461  for (unsigned j=0;j<num_nod;j++)
462  {
463  Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
464  double u;
465  Vector<double> x(2);
466  x[0]=nod_pt->x(0);
467  x[1]=nod_pt->x(1);
469  nod_pt->set_value(0,u);
470  }
471  }
472  }
473 } // end of actions_before_implicit_timestep
474 
475 
476 //=========start_of_actions_before_adapt==================================
477 /// Actions before adapt: Wipe the mesh of prescribed flux elements
478 //========================================================================
479 template<class ELEMENT>
481 {
482 
483  // Kill the flux elements and wipe surface mesh
484  delete_flux_elements(Surface_mesh_pt);
485 
486  // Rebuild the global mesh.
487  rebuild_global_mesh();
488 
489 } // end of actions_before_adapt
490 
491 
492 //==========start_of_actions_after_adapt==================================
493 /// Actions after adapt: Rebuild the mesh of prescribed flux elements
494 //========================================================================
495 template<class ELEMENT>
497 {
498  // Create prescribed-flux elements from all elements that are
499  // adjacent to boundary 0 and add them to surface mesh
500  create_flux_elements(0,Bulk_mesh_pt,Surface_mesh_pt);
501 
502  // Rebuild the global mesh
503  rebuild_global_mesh();
504 
505  // Loop over the flux elements to pass pointer to prescribed flux function
506  unsigned n_element=Surface_mesh_pt->nelement();
507  for(unsigned e=0;e<n_element;e++)
508  {
509  // Upcast from GeneralisedElement to UnsteadyHeat flux element
510  UnsteadyHeatFluxElement<ELEMENT> *el_pt =
511  dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*>(
512  Surface_mesh_pt->element_pt(e));
513 
514  // Set the pointer to the prescribed flux function
515  el_pt->flux_fct_pt() =
517  }
518 } // end of actions_after_adapt
519 
520 
521 //======================start_of_set_initial_condition====================
522 /// Set initial condition: Assign previous and current values
523 /// from exact solution.
524 //========================================================================
525 template<class ELEMENT>
527 {
528 
529  // Pointer to restart file
530  ifstream* restart_file_pt=0;
531 
532  // Restart?
533  //---------
534  // Restart file specified via command line [all programs have at least
535  // a single command line argument: their name. Ignore this here.]
536  if (CommandLineArgs::Argc==1)
537  {
538  cout << "No restart -- setting IC from exact solution" << std::endl;
539  }
540  else if (CommandLineArgs::Argc==2)
541  {
542  // Open restart file
543  restart_file_pt= new ifstream(CommandLineArgs::Argv[1],ios_base::in);
544  if (restart_file_pt!=0)
545  {
546  cout << "Have opened " << CommandLineArgs::Argv[1] <<
547  " for restart. " << std::endl;
548  }
549  else
550  {
551  std::ostringstream error_stream;
552  error_stream
553  << "ERROR while trying to open " << CommandLineArgs::Argv[1] <<
554  " for restart." << std::endl;
555 
556  throw OomphLibError(
557  error_stream.str(),
558  OOMPH_CURRENT_FUNCTION,
559  OOMPH_EXCEPTION_LOCATION);
560  }
561  }
562  // More than one command line argument?
563  else
564  {
565  std::ostringstream error_stream;
566  error_stream << "Can only specify one input file\n"
567  << "You specified the following command line arguments:\n";
568  //Fix this
569  CommandLineArgs::output();
570 
571  throw OomphLibError(
572  error_stream.str(),
573  OOMPH_CURRENT_FUNCTION,
574  OOMPH_EXCEPTION_LOCATION);
575  }
576 
577 
578  // Read restart data:
579  //-------------------
580  if (restart_file_pt!=0)
581  {
582  // Read the data from restart file and find out if the restart file
583  // was from an unsteady run
584  restart(*restart_file_pt);
585  }
586  // Assign initial condition from exact solution
587  //---------------------------------------------
588  else
589  {
590  // Choose initial timestep
591  double dt0=0.01;
592 
593  // Initialise timestep -- also sets the weights for all timesteppers
594  // in the problem.
595  initialise_dt(dt0);
596 
597  // Backup time in global timestepper
598  double backed_up_time=time_pt()->time();
599 
600  // Past history for velocities must be established for t=time0-deltat, ...
601  // Then provide current values (at t=time0) which will also form
602  // the initial guess for first solve at t=time0+deltat
603 
604  // Vector of exact solution value
605  Vector<double> soln(1);
606  Vector<double> x(2);
607 
608  //Find number of nodes in mesh
609  unsigned num_nod = Bulk_mesh_pt->nnode();
610 
611  // Get continuous times at previous timesteps
612  int nprev_steps=time_stepper_pt()->nprev_values();
613  Vector<double> prev_time(nprev_steps+1);
614  for (int itime=nprev_steps;itime>=0;itime--)
615  {
616  prev_time[itime]=time_pt()->time(unsigned(itime));
617  }
618 
619  // Loop over current & previous timesteps
620  for (int itime=nprev_steps;itime>=0;itime--)
621  {
622  double time=prev_time[itime];
623 
624  cout << "setting IC at time =" << time << std::endl;
625 
626  // Loop over the nodes to set initial guess everywhere
627  for (unsigned jnod=0;jnod<num_nod;jnod++)
628  {
629  // Get nodal coordinates
630  x[0]=Bulk_mesh_pt->node_pt(jnod)->x(0);
631  x[1]=Bulk_mesh_pt->node_pt(jnod)->x(1);
632 
633  // Get intial solution
635 
636  // Assign solution
637  Bulk_mesh_pt->node_pt(jnod)->set_value(itime,0,soln[0]);
638 
639  // Loop over coordinate directions
640  for (unsigned i=0;i<2;i++)
641  {
642  Bulk_mesh_pt->node_pt(jnod)->x(itime,i)=x[i];
643  }
644  }
645  }
646 
647  // Reset backed up time for global timestepper
648  time_pt()->time()=backed_up_time;
649  }
650 
651 } // end of set_initial_condition
652 
653 
654 //=======start_of_doc_solution============================================
655 /// Doc the solution
656 //========================================================================
657 template<class ELEMENT>
659 {
660  ofstream some_file;
661  char filename[100];
662 
663  // Number of plot points
664  unsigned npts;
665  npts=5;
666 
667  cout << std::endl;
668  cout << "=================================================" << std::endl;
669  cout << "Docing solution for t=" << time_pt()->time() << std::endl;
670  cout << "=================================================" << std::endl;
671 
672  // Output solution
673  //-----------------
674  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
675  Doc_info.number());
676  some_file.open(filename);
677  Bulk_mesh_pt->output(some_file,npts);
678  some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
679  << time_pt()->time() << "\"";
680  some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
681  some_file << "1" << std::endl;
682  some_file << "2" << std::endl;
683  some_file << " 0 0" << std::endl;
684  some_file << time_pt()->time()*20.0 << " 0" << std::endl;
685 
686  // Write dummy zones
687  some_file << "ZONE I=2,J=2" << std::endl;
688  some_file << "0.0 0.0 -1.2" << std::endl;
689  some_file << "1.3 0.0 -1.2" << std::endl;
690  some_file << "0.0 1.3 -1.2" << std::endl;
691  some_file << "1.3 1.3 -1.2" << std::endl;
692  some_file << "ZONE I=2,J=2" << std::endl;
693  some_file << "0.0 0.0 1.2" << std::endl;
694  some_file << "1.3 0.0 1.2" << std::endl;
695  some_file << "0.0 1.3 1.2" << std::endl;
696  some_file << "1.3 1.3 1.2" << std::endl;
697 
698  some_file.close();
699 
700 
701  // Output exact solution
702  //----------------------
703  sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
704  Doc_info.number());
705  some_file.open(filename);
706  Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
708 
709  // Write dummy zones
710  some_file << "ZONE I=2,J=2" << std::endl;
711  some_file << "0.0 0.0 -1.2" << std::endl;
712  some_file << "1.3 0.0 -1.2" << std::endl;
713  some_file << "0.0 1.3 -1.2" << std::endl;
714  some_file << "1.3 1.3 -1.2" << std::endl;
715  some_file << "ZONE I=2,J=2" << std::endl;
716  some_file << "0.0 0.0 1.2" << std::endl;
717  some_file << "1.3 0.0 1.2" << std::endl;
718  some_file << "0.0 1.3 1.2" << std::endl;
719  some_file << "1.3 1.3 1.2" << std::endl;
720 
721  some_file.close();
722 
723 
724  // Doc error
725  //----------
726  double error,norm;
727  sprintf(filename,"%s/error%i.dat",Doc_info.directory().c_str(),
728  Doc_info.number());
729  some_file.open(filename);
730  Bulk_mesh_pt->compute_error(some_file,
732  time_pt()->time(),
733  error,norm);
734  some_file.close();
735 
736 
737 
738  // Doc error and write trace file
739  //--------------------------------
740  cout << "error: " << error << std::endl;
741  cout << "norm : " << norm << std::endl << std::endl;
742 
743  Vector<double> x(2);
744  x[0]=Doc_node_pt->x(0);
745  x[1]=Doc_node_pt->x(1);
746  double u_exact;
747  TanhSolnForUnsteadyHeat::get_exact_u(time_pt()->time(),x,u_exact);
748  Vector<double > xi_wall(1);
749  Vector<double > r_wall(2);
750  xi_wall[0]=0.0;
751  Boundary_pt->position(xi_wall,r_wall);
752  Trace_file << time_pt()->time()
753  << " " << Doc_node_pt->value(0)
754  << " " << u_exact
755  << " " << r_wall[0]
756  << " " << TanhSolnForUnsteadyHeat::step_position(time_pt()->time())
757  << " " << Bulk_mesh_pt->nelement()
758  << " " << Bulk_mesh_pt->nrefined()
759  << " " << Bulk_mesh_pt->nunrefined()
760  << " " << error
761  << " " << norm << std::endl;
762 
763  // Plot wall posn
764  //---------------
765  sprintf(filename,"%s/Wall%i.dat",Doc_info.directory().c_str(),
766  Doc_info.number());
767  some_file.open(filename);
768 
769  unsigned nplot=100;
770  for (unsigned iplot=0;iplot<nplot;iplot++)
771  {
772  xi_wall[0]=0.5*Pi*double(iplot)/double(nplot-1);
773  Boundary_pt->position(xi_wall,r_wall);
774  some_file << r_wall[0] << " " << r_wall[1] << std::endl;
775  }
776  some_file.close();
777 
778  // Write restart file
779  sprintf(filename,"%s/restart%i.dat",Doc_info.directory().c_str(),
780  Doc_info.number());
781  some_file.open(filename);
782  dump_it(some_file);
783  some_file.close();
784 
785  // Increment number of doc
786  Doc_info.number()++;
787 
788 
789 } // end of doc_solution
790 
791 
792 //============start_of_create_flux_elements==============================
793 /// Create UnsteadyHeat Flux Elements on the b-th boundary of the Mesh object
794 /// pointed to by bulk_mesh_pt and add the elements to the Mesh object
795 /// pointed to by surface_mesh_pt.
796 //=======================================================================
797 template<class ELEMENT>
799 create_flux_elements(const unsigned& b, Mesh* const &bulk_mesh_pt,
800  Mesh* const &surface_mesh_pt)
801 {
802  // How many bulk elements are adjacent to boundary b?
803  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
804 
805  // Loop over the bulk elements adjacent to boundary b?
806  for(unsigned e=0;e<n_element;e++)
807  {
808  // Get pointer to the bulk element that is adjacent to boundary b
809  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
810  bulk_mesh_pt->boundary_element_pt(b,e));
811 
812  //Find the index of the face of element e along boundary b
813  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
814 
815  // Build the corresponding prescribed-flux element
816  UnsteadyHeatFluxElement<ELEMENT>* flux_element_pt = new
817  UnsteadyHeatFluxElement<ELEMENT>(bulk_elem_pt,face_index);
818 
819  //Add the prescribed-flux element to the surface mesh
820  surface_mesh_pt->add_element_pt(flux_element_pt);
821 
822  } //end of loop over bulk elements adjacent to boundary b
823 
824 } // end of create_flux_elements
825 
826 
827 //============start_of_delete_flux_elements==============================
828 /// Delete UnsteadyHeat Flux Elements and wipe the surface mesh
829 //=======================================================================
830 template<class ELEMENT>
832 delete_flux_elements(Mesh* const &surface_mesh_pt)
833 {
834  // How many surface elements are in the surface mesh
835  unsigned n_element = surface_mesh_pt->nelement();
836 
837  // Loop over the surface elements
838  for(unsigned e=0;e<n_element;e++)
839  {
840  // Kill surface element
841  delete surface_mesh_pt->element_pt(e);
842  }
843 
844  // Wipe the mesh
845  surface_mesh_pt->flush_element_and_node_storage();
846 
847 } // end of delete_flux_elements
848 
849 
850 //=======start_of_dump_it=================================================
851 /// Dump the solution to disk
852 //========================================================================
853 template<class ELEMENT>
855 {
856  // Dump the refinement pattern and the generic problem data
857  Problem::dump(dump_file);
858 
859 } // end of dump_it
860 
861 //=========start_of_restart===============================================
862 /// Read solution from disk
863 //========================================================================
864 template<class ELEMENT>
866 {
867 
868  // Refine the mesh and read in the generic problem data
869  Problem::read(restart_file);
870 
871 } // end of restart
872 
873 
874 
875 /// /////////////////////////////////////////////////////////////////////
876 /// /////////////////////////////////////////////////////////////////////
877 /// /////////////////////////////////////////////////////////////////////
878 
879 //======start_of_main=====================================================
880 /// Demonstrate how to solve an unsteady heat problem
881 /// with mesh adaptation. Command line arguments specify
882 /// the name of the restart file.
883 //========================================================================
884 int main(int argc, char* argv[])
885 {
886 
887  // Store command line arguments
888  CommandLineArgs::setup(argc,argv);
889 
890  // Build problem
893 
894  // Specify duration of the simulation
895  //double t_max=3.0;
896 
897  // Set targets for spatial adaptivity
898  problem.bulk_mesh_pt()->max_permitted_error()=0.001;
899  problem.bulk_mesh_pt()->min_permitted_error()=0.0001;
900 
901  // First timestep?
902  bool first=true;
903 
904  // Max. number of spatial adaptations per timestep. Allow plenty
905  // of adaptations at first timestep as the initial conditions
906  // can be reset "exactly" from without any interpolation error.
907  unsigned max_adapt=10;
908 
909  // Set IC
910  problem.set_initial_condition();
911 
912  // Initial timestep: Use the one used when setting up the initial
913  // condition
914  double dt=problem.time_pt()->dt();
915 
916  // If restart: The first step isn't really the first step,
917  // i.e. initial condition should not be re-set when
918  // adaptive refinement has been performed. Also, limit
919  // the max. number of refinements per timestep to the
920  // normal value straightaway.
921  if (CommandLineArgs::Argc==2)
922  {
923  first=false;
924  max_adapt=1;
925  }
926  // If no restart, refine mesh uniformly before we get started
927  else
928  {
929  problem.refine_uniformly();
930  problem.refine_uniformly();
931  // Solution is automatically interpolated from the coarse initial mesh
932  // onto the refined mesh but this provides a very poor representation
933  // of the intial condition: Re-assign the initial conditions
934  problem.set_initial_condition();
935  }
936 
937  //Output FE representation of the initial condition
938  problem.doc_solution();
939 
940  // Find number of steps
941  unsigned nstep = 6; // unsigned(t_max/dt);
942 
943  // Timestepping loop
944  for (unsigned istep=0;istep<nstep;istep++)
945  {
946  // Take timestep
947  problem.unsteady_newton_solve(dt,max_adapt,first);
948 
949  // Now we've done the first timestep -- don't re-set the IC
950  // in subsequent steps
951  first=false;
952 
953  // Reduce the number of spatial adaptations to one per
954  // timestep
955  max_adapt=1;
956 
957  //Output solution
958  problem.doc_solution();
959 
960  }
961 
962 
963 }; // end of main
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
MyUnitCircle()
Constructor: The circle is a 1D object (i.e. it's parametrised by one intrinsic coordinate) in 2D spa...
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Parametrised position on object: r(xi). Evaluated at previous time level. t=0: current time; t>0: pre...
void position(const Vector< double > &xi, Vector< double > &r) const
Current position vector to material point at Lagrangian coordinate xi.
virtual ~MyUnitCircle()
Destructor: Empty.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed flux elements.
void dump_it(ofstream &dump_file)
Dump problem data to allow for later restart.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function. Note that his overlo...
GeomObject * Boundary_pt
Pointer to GeomObject that specifies the domain bondary.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed flux elements.
RefineableQuarterCircleSectorMesh< ELEMENT > * bulk_mesh_pt()
Pointer to bulk mesh.
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &surface_mesh_pt)
Create UnsteadyHeat flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them t...
Node * Doc_node_pt
Pointer to central node (exists at all refinement levels) for doc.
RefineableUnsteadyHeatProblem(UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
~RefineableUnsteadyHeatProblem()
Destructor: Close trace file.
UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt Source_fct_pt
Pointer to source function.
void actions_after_implicit_timestep()
Update the problem specs after timestep (empty)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Set Dirchlet boundary conditions from exact solution.
void delete_flux_elements(Mesh *const &surface_mesh_pt)
Delete UnsteadyHeat flux elements and wipe the surface mesh.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void restart(ifstream &restart_file)
Read problem data for restart.
RefineableQuarterCircleSectorMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void get_exact_u(const double &time, const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double Alpha
Parameter for steepness of step.
double Gamma
Parameter for timescale of step translation.
void get_exact_u(const double &time, const Vector< double > &x, double &u)
Exact solution as a scalar.
double Beta
Parameter for amplitude of step translation.
double step_position(const double &time)
Position of step (x-axis intercept)
void get_source(const double &time, const Vector< double > &x, double &source)
Source function to make it an exact solution.
void prescribed_flux_on_fixed_y_boundary(const double &time, const Vector< double > &x, double &flux)
Flux required by the exact solution on a boundary on which y is fixed.
double TanPhi
Parameter for angle of step.
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...