32 #include "unsteady_heat.h"
35 #include "meshes/quarter_circle_sector_mesh.h"
39 using namespace oomph;
41 using namespace MathematicalConstants;
68 void position(
const Vector<double>& xi, Vector<double>& r)
const
79 void position(
const unsigned& t,
const Vector<double>& xi,
80 Vector<double>& r)
const
117 return Beta*tanh(
Gamma*cos(0.2E1*MathematicalConstants::Pi*time));
127 MathematicalConstants::Pi*time)))-Y));
131 void get_exact_u(
const double& time,
const Vector<double>& x,
double& u)
136 MathematicalConstants::Pi*time)))-Y));
140 void get_source(
const double& time,
const Vector<double>& x,
double& source)
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*
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*
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;
158 const Vector<double>& x,
170 MathematicalConstants::Pi*time)))-Y)),2.0))*
Alpha*
TanPhi*NX+(1.0-pow(tanh(
172 time)))-Y)),2.0))*
Alpha*NY;
186 template<
class ELEMENT>
194 UnsteadyHeatSourceFctPt source_fct_pt);
210 void actions_before_implicit_timestep();
213 void actions_before_adapt();
216 void actions_after_adapt();
223 void set_initial_condition();
228 void create_flux_elements(
const unsigned &b, Mesh*
const &bulk_mesh_pt,
229 Mesh*
const &surface_mesh_pt);
232 void delete_flux_elements(Mesh*
const &surface_mesh_pt);
238 void dump_it(ofstream& dump_file);
241 void restart(ifstream& restart_file);
279 template<
class ELEMENT>
281 UnsteadyHeatEquations<2>::UnsteadyHeatSourceFctPt source_fct_pt) :
282 Source_fct_pt(source_fct_pt)
297 sprintf(filename,
"%s/trace.dat",
Doc_info.directory().c_str());
300 Trace_file <<
"VARIABLES=\"time t\",\"u<SUB>FE</SUB>\",\"u<SUB>exact</SUB>\","
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\""
328 add_time_stepper_pt(
new BDF<2>());
341 double xi_hi=MathematicalConstants::Pi/2.0;
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());
366 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
367 Bulk_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
373 for(
unsigned b=0;b<n_bound;b++)
380 for (
unsigned n=0;n<n_node;n++)
389 FiniteElement* el0_pt=
Bulk_mesh_pt->finite_element_pt(0);
390 unsigned nnod=el0_pt->nnode();
402 for(
unsigned i=0;i<n_element;i++)
405 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(
Bulk_mesh_pt->element_pt(i));
413 for(
unsigned e=0;e<n_element;e++)
416 UnsteadyHeatFluxElement<ELEMENT> *el_pt =
417 dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*
>(
421 el_pt->flux_fct_pt() =
426 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
434 template<
class ELEMENT>
446 template<
class ELEMENT>
450 double time=time_pt()->time();
453 unsigned num_bound = Bulk_mesh_pt->nboundary();
454 for(
unsigned b=0;b<num_bound;b++)
460 unsigned num_nod=Bulk_mesh_pt->nboundary_node(b);
461 for (
unsigned j=0;j<num_nod;j++)
463 Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(b,j);
469 nod_pt->set_value(0,u);
479 template<
class ELEMENT>
484 delete_flux_elements(Surface_mesh_pt);
487 rebuild_global_mesh();
495 template<
class ELEMENT>
500 create_flux_elements(0,Bulk_mesh_pt,Surface_mesh_pt);
503 rebuild_global_mesh();
506 unsigned n_element=Surface_mesh_pt->nelement();
507 for(
unsigned e=0;e<n_element;e++)
510 UnsteadyHeatFluxElement<ELEMENT> *el_pt =
511 dynamic_cast<UnsteadyHeatFluxElement<ELEMENT>*
>(
512 Surface_mesh_pt->element_pt(e));
515 el_pt->flux_fct_pt() =
525 template<
class ELEMENT>
530 ifstream* restart_file_pt=0;
536 if (CommandLineArgs::Argc==1)
538 cout <<
"No restart -- setting IC from exact solution" << std::endl;
540 else if (CommandLineArgs::Argc==2)
543 restart_file_pt=
new ifstream(CommandLineArgs::Argv[1],ios_base::in);
544 if (restart_file_pt!=0)
546 cout <<
"Have opened " << CommandLineArgs::Argv[1] <<
547 " for restart. " << std::endl;
551 std::ostringstream error_stream;
553 <<
"ERROR while trying to open " << CommandLineArgs::Argv[1] <<
554 " for restart." << std::endl;
558 OOMPH_CURRENT_FUNCTION,
559 OOMPH_EXCEPTION_LOCATION);
565 std::ostringstream error_stream;
566 error_stream <<
"Can only specify one input file\n"
567 <<
"You specified the following command line arguments:\n";
569 CommandLineArgs::output();
573 OOMPH_CURRENT_FUNCTION,
574 OOMPH_EXCEPTION_LOCATION);
580 if (restart_file_pt!=0)
584 restart(*restart_file_pt);
598 double backed_up_time=time_pt()->time();
605 Vector<double> soln(1);
609 unsigned num_nod = Bulk_mesh_pt->nnode();
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--)
616 prev_time[itime]=time_pt()->time(
unsigned(itime));
620 for (
int itime=nprev_steps;itime>=0;itime--)
622 double time=prev_time[itime];
624 cout <<
"setting IC at time =" << time << std::endl;
627 for (
unsigned jnod=0;jnod<num_nod;jnod++)
630 x[0]=Bulk_mesh_pt->node_pt(jnod)->x(0);
631 x[1]=Bulk_mesh_pt->node_pt(jnod)->x(1);
637 Bulk_mesh_pt->node_pt(jnod)->set_value(itime,0,soln[0]);
640 for (
unsigned i=0;i<2;i++)
642 Bulk_mesh_pt->node_pt(jnod)->x(itime,i)=x[i];
648 time_pt()->time()=backed_up_time;
657 template<
class ELEMENT>
668 cout <<
"=================================================" << std::endl;
669 cout <<
"Docing solution for t=" << time_pt()->time() << std::endl;
670 cout <<
"=================================================" << std::endl;
674 sprintf(filename,
"%s/soln%i.dat",Doc_info.directory().c_str(),
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;
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;
703 sprintf(filename,
"%s/exact_soln%i.dat",Doc_info.directory().c_str(),
705 some_file.open(filename);
706 Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
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;
727 sprintf(filename,
"%s/error%i.dat",Doc_info.directory().c_str(),
729 some_file.open(filename);
730 Bulk_mesh_pt->compute_error(some_file,
740 cout <<
"error: " << error << std::endl;
741 cout <<
"norm : " << norm << std::endl << std::endl;
744 x[0]=Doc_node_pt->x(0);
745 x[1]=Doc_node_pt->x(1);
748 Vector<double > xi_wall(1);
749 Vector<double > r_wall(2);
751 Boundary_pt->position(xi_wall,r_wall);
752 Trace_file << time_pt()->time()
753 <<
" " << Doc_node_pt->value(0)
757 <<
" " << Bulk_mesh_pt->nelement()
758 <<
" " << Bulk_mesh_pt->nrefined()
759 <<
" " << Bulk_mesh_pt->nunrefined()
761 <<
" " << norm << std::endl;
765 sprintf(filename,
"%s/Wall%i.dat",Doc_info.directory().c_str(),
767 some_file.open(filename);
770 for (
unsigned iplot=0;iplot<nplot;iplot++)
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;
779 sprintf(filename,
"%s/restart%i.dat",Doc_info.directory().c_str(),
781 some_file.open(filename);
797 template<
class ELEMENT>
800 Mesh*
const &surface_mesh_pt)
803 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
806 for(
unsigned e=0;e<n_element;e++)
809 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
810 bulk_mesh_pt->boundary_element_pt(b,e));
813 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
816 UnsteadyHeatFluxElement<ELEMENT>* flux_element_pt =
new
817 UnsteadyHeatFluxElement<ELEMENT>(bulk_elem_pt,face_index);
820 surface_mesh_pt->add_element_pt(flux_element_pt);
830 template<
class ELEMENT>
835 unsigned n_element = surface_mesh_pt->nelement();
838 for(
unsigned e=0;e<n_element;e++)
841 delete surface_mesh_pt->element_pt(e);
845 surface_mesh_pt->flush_element_and_node_storage();
853 template<
class ELEMENT>
857 Problem::dump(dump_file);
864 template<
class ELEMENT>
869 Problem::read(restart_file);
884 int main(
int argc,
char* argv[])
888 CommandLineArgs::setup(argc,argv);
907 unsigned max_adapt=10;
914 double dt=problem.time_pt()->dt();
921 if (CommandLineArgs::Argc==2)
929 problem.refine_uniformly();
930 problem.refine_uniformly();
944 for (
unsigned istep=0;istep<nstep;istep++)
947 problem.unsteady_newton_solve(dt,max_adapt,first);
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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...
void doc_solution()
Doc the solution.
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.
ofstream Trace_file
Trace file.
~RefineableUnsteadyHeatProblem()
Destructor: Close trace file.
UnsteadyHeatEquations< 2 >::UnsteadyHeatSourceFctPt Source_fct_pt
Pointer to source function.
DocInfo Doc_info
Doc info object.
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[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...