35 #include "meshes/simple_rectangular_quadmesh.h"
39 using namespace oomph;
73 double N[2] = {1.0, 0.0};
92 template<
class ELEMENT>
106 void doc_solution(DocInfo& doc_info);
112 void actions_before_newton_solve();
119 void create_flux_elements(
const unsigned &b);
136 template<
class ELEMENT>
139 : Source_fct_pt(source_fct_pt)
157 Problem::mesh_pt()=
new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
169 unsigned n_bound = mesh_pt()->nboundary();
170 for(
unsigned b=0;b<n_bound;b++)
175 unsigned n_node= mesh_pt()->nboundary_node(b);
176 for (
unsigned n=0;n<n_node;n++)
178 mesh_pt()->boundary_node_pt(b,n)->pin(0);
189 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(e));
196 unsigned n_element=mesh_pt()->nelement();
203 PoissonFluxElement<ELEMENT> *el_pt =
204 dynamic_cast<PoissonFluxElement<ELEMENT>*
>(mesh_pt()->element_pt(e));
207 el_pt->flux_fct_pt() =
212 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
220 template<
class ELEMENT>
224 unsigned n_element = mesh_pt()->nboundary_element(b);
227 for(
unsigned e=0;e<n_element;e++)
230 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
231 mesh_pt()->boundary_element_pt(b,e));
234 int face_index = mesh_pt()->face_index_at_boundary(b,e);
237 PoissonFluxElement<ELEMENT>* flux_element_pt =
new
238 PoissonFluxElement<ELEMENT>(bulk_elem_pt,face_index);
241 mesh_pt()->add_element_pt(flux_element_pt);
253 template<
class ELEMENT>
257 unsigned n_bound = mesh_pt()->nboundary();
260 for(
unsigned i=0;i<n_bound;i++)
266 unsigned n_node = mesh_pt()->nboundary_node(i);
269 for (
unsigned n=0;n<n_node;n++)
272 Node* nod_pt = mesh_pt()->boundary_node_pt(i,n);
284 nod_pt->set_value(0,u[0]);
295 template<
class ELEMENT>
308 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
310 some_file.open(filename);
311 mesh_pt()->output(some_file,npts);
316 sprintf(filename,
"%s/exact_soln%i.dat",doc_info.directory().c_str(),
318 some_file.open(filename);
319 for(
unsigned e=0;e<Npoisson_elements;e++)
321 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(e));
356 doc_info.set_directory(
"RESLT");
365 cout <<
"\n\n\nProblem self-test ";
366 if (problem.self_test()==0)
368 cout <<
"passed: Problem can be solved." << std::endl;
372 throw OomphLibError(
"Self test failed",
373 OOMPH_CURRENT_FUNCTION,
374 OOMPH_EXCEPTION_LOCATION);
387 for (
unsigned istep=0;istep<nstep;istep++)
392 cout <<
"\n\nSolving for TanhSolnForPoisson::Alpha="
396 problem.newton_solve();
2D Poisson problem on rectangular domain, discretised with 2D QPoisson elements. Flux boundary condit...
FluxPoissonProblem(PoissonEquations< 2 >::PoissonSourceFctPt source_fct_pt)
Constructor: Pass pointer to source function.
unsigned Npoisson_elements
Number of Poisson "bulk" elements (We're attaching the flux elements to the bulk mesh --> only the fi...
void create_flux_elements(const unsigned &b)
Create Poisson flux elements on the b-th boundary of the problem's mesh.
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Reset boundary conditions to the values from the exact solutio...
PoissonEquations< 2 >::PoissonSourceFctPt Source_fct_pt
Pointer to source function.
void doc_solution(DocInfo &doc_info)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to.
~FluxPoissonProblem()
Destructor (empty)
Namespace for exact solution for Poisson equation with "sharp step".
void prescribed_flux_on_fixed_x_boundary(const Vector< double > &x, double &flux)
Flux required by the exact solution on a boundary on which x is fixed.
double TanPhi
Parameter for angle Phi of "step".
void source_function(const Vector< double > &x, double &source)
Source function required to make the solution above an exact solution.
double Alpha
Parameter for steepness of "step".
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
int main()
Demonstrate how to solve 2D Poisson problem with flux boundary conditions.