32 #include "advection_diffusion.h"
35 #include "meshes/rectangular_quadmesh.h"
39 using namespace oomph;
75 wind[0]=sin(6.0*x[1]);
76 wind[1]=cos(6.0*x[0]);
92 template<
class ELEMENT>
101 AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
102 AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt,
103 const bool& use_stabilisation);
111 void actions_before_newton_solve();
124 return dynamic_cast<RectangularQuadMesh<ELEMENT>*
>(
134 AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt
Source_fct_pt;
137 AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt
Wind_fct_pt;
151 template<
class ELEMENT>
153 AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,
154 AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt,
155 const bool& use_stabilisation)
156 : Source_fct_pt(source_fct_pt), Wind_fct_pt(wind_fct_pt),
157 Use_stabilisation(use_stabilisation)
161 if (use_stabilisation)
163 Doc_info.set_directory(
"RESLT_stabilised");
167 Doc_info.set_directory(
"RESLT_unstabilised");
186 new RectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
192 unsigned num_bound =
mesh_pt()->nboundary();
193 for(
unsigned ibound=0;ibound<num_bound;ibound++)
195 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
196 for (
unsigned inod=0;inod<num_nod;inod++)
198 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
207 unsigned n_element =
mesh_pt()->nelement();
208 for(
unsigned i=0;i<n_element;i++)
211 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(
mesh_pt()->element_pt(i));
224 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
232 template<
class ELEMENT>
236 unsigned num_bound = mesh_pt()->nboundary();
239 for(
unsigned ibound=0;ibound<num_bound;ibound++)
242 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
245 for (
unsigned inod=0;inod<num_nod;inod++)
248 Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);
260 nod_pt->set_value(0,u[0]);
265 unsigned n_element = mesh_pt()->nelement();
266 for(
unsigned i=0;i<n_element;i++)
269 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
273 if (Use_stabilisation)
276 el_pt->compute_stabilisation_parameter();
281 el_pt->switch_off_stabilisation();
293 template<
class ELEMENT>
305 sprintf(filename,
"%s/soln%i.dat",Doc_info.directory().c_str(),
307 some_file.open(filename);
308 mesh_pt()->output(some_file,npts);
322 bool use_stabilisation=
true;
333 problem.newton_solve();
345 bool use_stabilisation=
false;
356 problem.newton_solve();
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
~SUPGAdvectionDiffusionProblem()
Destructor. Empty.
void doc_solution()
Doc the solution.
void actions_after_newton_solve()
Update the problem after solve (empty)
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function.
SUPGAdvectionDiffusionProblem(AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt source_fct_pt, AdvectionDiffusionEquations< 2 >::AdvectionDiffusionWindFctPt wind_fct_pt, const bool &use_stabilisation)
Constructor: Pass pointer to source and wind functions, and flag to indicate if stabilisation is to b...
AdvectionDiffusionEquations< 2 >::AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function.
bool Use_stabilisation
Flag to indicate if stabilisation is to be used.
RectangularQuadMesh< ELEMENT > * mesh_pt()
Overloaded version of the problem's access function to the mesh. Recasts the pointer to the base Mesh...
DocInfo Doc_info
DocInfo object.
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the tanh solution...
Namespace for global parameters: Unforced problem with boundary values corresponding to a steep tanh ...
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
void get_boundary_values(const Vector< double > &x, Vector< double > &u)
Some "solution" for assignment of boundary values.
void source_function(const Vector< double > &x_vect, double &source)
Zero source function.
double Peclet
Peclet number.
double Alpha
Parameter for steepness of step in boundary values.
double TanPhi
Parameter for angle of step in boundary values: 45 degrees.
int main()
Driver code for 2D AdvectionDiffusion problem.