33 #include "linear_wave.h"
36 #include "meshes/rectangular_quadmesh.h"
40 using namespace oomph;
56 double exact_u(
const double& time,
const Vector<double>& x)
58 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
59 return tanh(1.0-
Alpha*(zeta-time));
63 double exact_dudt(
const double& time,
const Vector<double>& x)
65 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
67 cosh(1.0-
Alpha*(zeta-time)));
71 double exact_d2udt2(
const double& time,
const Vector<double>& x)
73 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
75 (cosh(1.0-
Alpha*(zeta-time))*cosh(1.0-
Alpha*(zeta-time)));
80 void get_exact_u(
const double& time,
const Vector<double>& x,
89 void get_source(
const double& time,
const Vector<double>& x,
double& source)
99 double zeta=cos(
Phi)*x[0]+sin(
Phi)*x[1];
101 *cosh(1.0-
Alpha*(zeta-time)));
103 *cosh(1.0-
Alpha*(zeta-time)));
109 const Vector<double>& x,
113 Vector<double> dudx(2), n(2);
121 flux = dudx[0]*n[0] + dudx[1]*n[1];
139 template<
class ELEMENT,
class TIMESTEPPER>
148 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt);
166 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
167 initial_value_fct(1);
168 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
169 initial_veloc_fct(1);
170 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
171 initial_accel_fct(1);
179 unsigned num_bound=Bulk_mesh_pt->nboundary();
180 for (
unsigned ibound=0;ibound<num_bound;ibound++)
186 unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
187 for (
unsigned inod=0;inod<num_nod;inod++)
190 Node* nod_pt=Bulk_mesh_pt->boundary_node_pt(ibound,inod);
192 bool use_direct_assignment=
false;
193 if (use_direct_assignment)
208 TIMESTEPPER* timestepper_pt=
dynamic_cast<TIMESTEPPER*
>
212 timestepper_pt->assign_initial_data_values(nod_pt,
237 void create_flux_elements(
const unsigned& b, Mesh*
const& bulk_mesh_pt,
238 Mesh*
const& surface_mesh_pt);
256 template<
class ELEMENT,
class TIMESTEPPER>
258 const unsigned& nx,
const unsigned& ny,
259 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt)
263 add_time_stepper_pt(
new TIMESTEPPER());
293 Bulk_mesh_pt=
new RectangularQuadMesh<ELEMENT>(Nx,Ny,Lx,Ly,time_stepper_pt());
298 Surface_mesh_pt =
new Mesh;
304 create_flux_elements(1,Bulk_mesh_pt,Surface_mesh_pt);
307 add_sub_mesh(Bulk_mesh_pt);
308 add_sub_mesh(Surface_mesh_pt);
318 unsigned num_bound = Bulk_mesh_pt->nboundary();
319 for(
unsigned ibound=0;ibound<num_bound;ibound++)
323 unsigned num_nod= Bulk_mesh_pt->nboundary_node(ibound);
324 for (
unsigned inod=0;inod<num_nod;inod++)
326 Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(0);
336 unsigned n_element = Bulk_mesh_pt->nelement();
337 for(
unsigned i=0;i<n_element;i++)
340 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(i));
343 el_pt->source_fct_pt() = source_fct_pt;
352 n_element = Surface_mesh_pt->nelement();
353 for (
unsigned e=0;e<n_element;e++)
356 LinearWaveFluxElement<ELEMENT> *el_pt =
357 dynamic_cast< LinearWaveFluxElement<ELEMENT>*
>(
358 Surface_mesh_pt->element_pt(e));
360 el_pt->flux_fct_pt() =
365 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
374 template<
class ELEMENT,
class TIMESTEPPER>
379 TIMESTEPPER* timestepper_pt=
dynamic_cast<TIMESTEPPER*
>(time_stepper_pt());
384 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
385 initial_value_fct(1);
386 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
387 initial_veloc_fct(1);
388 Vector<typename TIMESTEPPER::NodeInitialConditionFctPt>
389 initial_accel_fct(1);
400 unsigned num_nod=mesh_pt()->nnode();
401 for (
unsigned jnod=0;jnod<num_nod;jnod++)
404 Node* nod_pt=mesh_pt()->node_pt(jnod);
407 timestepper_pt->assign_initial_data_values(nod_pt,
421 template<
class ELEMENT,
class TIMESTEPPER>
424 Mesh*
const &surface_mesh_pt)
427 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
430 for(
unsigned e=0;e<n_element;e++)
433 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
434 bulk_mesh_pt->boundary_element_pt(b,e));
437 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
440 LinearWaveFluxElement<ELEMENT>* flux_element_pt =
new
441 LinearWaveFluxElement<ELEMENT>(bulk_elem_pt,face_index);
444 surface_mesh_pt->add_element_pt(flux_element_pt);
455 template<
class ELEMENT,
class TIMESTEPPER>
467 cout <<
"=================================================" << std::endl;
468 cout <<
"Docing solution for t=" << time_pt()->time() << std::endl;
469 cout <<
"=================================================" << std::endl;
473 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
475 some_file.open(filename);
476 Bulk_mesh_pt->output(some_file,npts);
477 some_file <<
"TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
478 << time_pt()->time() <<
"\"";
479 some_file <<
"GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
480 some_file <<
"1" << std::endl;
481 some_file <<
"2" << std::endl;
482 some_file <<
" 0 0" << std::endl;
483 some_file << time_pt()->time()*20.0 <<
" 0" << std::endl;
488 sprintf(filename,
"%s/exact_soln%i.dat",doc_info.directory().c_str(),
490 some_file.open(filename);
491 Bulk_mesh_pt->output_fct(some_file,npts,time_pt()->time(),
498 sprintf(filename,
"%s/error%i.dat",doc_info.directory().c_str(),
500 some_file.open(filename);
501 Bulk_mesh_pt->compute_error(some_file,
507 cout <<
"error: " << error << std::endl;
508 cout <<
"norm : " << norm << std::endl << std::endl;
509 Trace_file << time_pt()->time() <<
" " << time_pt()->dt()
510 <<
" " << Bulk_mesh_pt->nelement() <<
" "
511 << error <<
" " << norm << std::endl;
521 template<
class ELEMENT,
class TIMESTEPPER>
530 doc_info.set_directory(
"RESLT");
537 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
538 Trace_file.open(filename);
549 time_pt()->time()=time0;
552 time_pt()->initialise_dt(dt);
555 set_initial_condition();
558 doc_solution(doc_info);
564 unsigned nstep=unsigned(t_max/dt);
567 if (CommandLineArgs::Argc>1)
570 cout <<
"validation run" << std::endl;
574 for (
unsigned istep=0;istep<nstep;istep++)
577 unsteady_newton_solve(dt);
580 doc_solution(doc_info);
603 int main(
int argc,
char* argv[])
607 CommandLineArgs::setup(argc,argv);
610 LinearWaveEquations<2>::LinearWaveSourceFctPt source_fct_pt=
621 problem(n_x,n_y,source_fct_pt);
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &surface_mesh_pt)
Create LinearWave flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to ...
void actions_before_implicit_timestep()
Update the problem specs before next timestep: Set Dirchlet boundary conditions from exact solution.
LinearWaveProblem(const unsigned &nx, const unsigned &ny, const bool &impulsive_start, LinearWaveEquations< 2 >::LinearWaveSourceFctPt source_fct_pt)
Constructor: pass number of elements in x and y directions, bool indicating impulsive or "smooth" sta...
void actions_after_implicit_timestep()
Update the problem specs after solve (empty)
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void doc_solution(DocInfo &doc_info)
Doc the solution.
void unsteady_run()
Do unsteady run.
RectangularQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
~LinearWaveProblem()
Destructor (empty)
void set_initial_condition()
Set initial condition (incl history values) according to specified function.
Namespace for exact solution for LinearWave equation with sharp step.
void get_source(const double &time, const Vector< double > &x, double &source)
Source function to make it an exact solution.
double Alpha
Parameter for steepness of step.
double Phi
Orientation of step wave.
void prescribed_flux_on_fixed_y_boundary(const double &time, const Vector< double > &x, double &flux)
Prescribed flux on a fixed y max boundary.
double exact_d2udt2(const double &time, const Vector< double > &x)
2nd time-deriv of exact solution
void get_exact_u(const double &time, const Vector< double > &x, Vector< double > &u)
Exact solution as a vector.
double exact_dudt(const double &time, const Vector< double > &x)
1st time-deriv of exact solution
double exact_u(const double &time, const Vector< double > &x)
Exact solution.
void get_exact_gradient(const double &time, const Vector< double > &x, Vector< double > &dudx)
Gradient of exact solution.
int main(int argc, char *argv[])
///////////////////////////////////////////////////////////////////// ///////////////////////////////...