35 #include "meshes/triangle_mesh.h"
38 using namespace oomph;
55 void get_exact_u(
const Vector<double>& x, Vector<double>& u)
61 void get_source(
const Vector<double>& x,
double& source)
72 void zero(
const Vector<double>& x, Vector<double>& u)
89 template<
class ELEMENT>
107 apply_boundary_conditions();
141 template<
class ELEMENT>
145 Vector<double> zeta(1);
148 Vector<double> posn(2);
151 double x_center = 0.0;
152 double y_center = 0.0;
155 Ellipse * outer_boundary_ellipse_pt =
new Ellipse(A,B);
158 TriangleMeshClosedCurve* closed_curve_pt=0;
162 bool polygon_for_outer_boundary=
false;
164 polygon_for_outer_boundary=
true;
166 if (polygon_for_outer_boundary)
170 double unit_zeta = 0.5*MathematicalConstants::Pi/double(n_seg);
174 Vector<TriangleMeshCurveSection*> boundary_polyline_pt(2);
177 Vector<Vector<double> > bound_coords(n_seg+1);
181 for(
unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
184 bound_coords[ipoint].resize(2);
187 zeta[0]=unit_zeta*double(ipoint);
188 outer_boundary_ellipse_pt->position(zeta,posn);
189 bound_coords[ipoint][0]=posn[0]+x_center;
190 bound_coords[ipoint][1]=posn[1]+y_center;
194 unsigned boundary_id=0;
195 boundary_polyline_pt[0]=
new TriangleMeshPolyLine(bound_coords,boundary_id);
202 for(
unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
205 bound_coords[ipoint].resize(2);
208 zeta[0]=(unit_zeta*double(ipoint))+0.5*MathematicalConstants::Pi;
209 outer_boundary_ellipse_pt->position(zeta,posn);
210 bound_coords[ipoint][0]=posn[0]+x_center;
211 bound_coords[ipoint][1]=posn[1]+y_center;
216 boundary_polyline_pt[1]=
new TriangleMeshPolyLine(bound_coords,boundary_id);
221 TriangleMeshPolygon* outer_boundary_polygon_pt =
222 new TriangleMeshPolygon(boundary_polyline_pt);
225 outer_boundary_polygon_pt->
226 enable_redistribution_of_segments_between_polylines();
229 closed_curve_pt=outer_boundary_polygon_pt;
238 Vector<TriangleMeshCurveSection*> outer_curvilinear_boundary_pt(2);
242 double zeta_start=0.0;
243 double zeta_end=MathematicalConstants::Pi;
245 unsigned boundary_id=0;
246 outer_curvilinear_boundary_pt[0]=
new TriangleMeshCurviLine(
247 outer_boundary_ellipse_pt,zeta_start,zeta_end,nsegment,boundary_id);
251 zeta_start=MathematicalConstants::Pi;
252 zeta_end=2.0*MathematicalConstants::Pi;
255 outer_curvilinear_boundary_pt[1]=
new TriangleMeshCurviLine(
256 outer_boundary_ellipse_pt,zeta_start,zeta_end,nsegment,boundary_id);
260 TriangleMeshClosedCurve* curvilinear_outer_boundary_pt=
261 new TriangleMeshClosedCurve(outer_curvilinear_boundary_pt);
264 closed_curve_pt=curvilinear_outer_boundary_pt;
270 Vector<TriangleMeshClosedCurve*> hole_pt(2);
280 Ellipse* polygon_ellipse_pt=
new Ellipse(A,B);
284 double unit_zeta = MathematicalConstants::Pi/double(n_seg);
288 Vector<TriangleMeshCurveSection*> hole_polyline_pt(2);
295 Vector<Vector<double> > bound_hole(n_seg+1);
296 for(
unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
299 bound_hole[ipoint].resize(2);
302 zeta[0]=unit_zeta*double(ipoint);
303 polygon_ellipse_pt->position(zeta,posn);
304 bound_hole[ipoint][0]=posn[0]+x_center;
305 bound_hole[ipoint][1]=posn[1]+y_center;
309 unsigned boundary_id=2;
312 hole_polyline_pt[0] =
new TriangleMeshPolyLine(bound_hole,boundary_id);
317 for(
unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
320 bound_hole[ipoint].resize(2);
323 zeta[0]=(unit_zeta*double(ipoint))+MathematicalConstants::Pi;
324 polygon_ellipse_pt->position(zeta,posn);
325 bound_hole[ipoint][0]=posn[0]+x_center;
326 bound_hole[ipoint][1]=posn[1]+y_center;
333 hole_polyline_pt[1] =
new TriangleMeshPolyLine(bound_hole,boundary_id);
340 Vector<double> hole_center(2);
341 hole_center[0]=x_center;
342 hole_center[1]=y_center;
344 hole_pt[0] =
new TriangleMeshPolygon(hole_polyline_pt, hole_center);
355 Ellipse* ellipse_pt=
new Ellipse(A,B);
358 Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt(2);
363 double zeta_start=0.0;
364 double zeta_end=MathematicalConstants::Pi;
365 unsigned nsegment=10;
367 curvilinear_boundary_pt[0]=
new TriangleMeshCurviLine(
368 ellipse_pt,zeta_start,zeta_end,
369 nsegment,boundary_id);
373 zeta_start=MathematicalConstants::Pi;
374 zeta_end=2.0*MathematicalConstants::Pi;
377 curvilinear_boundary_pt[1]=
new TriangleMeshCurviLine(
378 ellipse_pt,zeta_start,zeta_end,
379 nsegment,boundary_id);
384 Vector<double> hole_coords(2);
387 Vector<TriangleMeshClosedCurve*> curvilinear_hole_pt(1);
389 new TriangleMeshClosedCurve(curvilinear_boundary_pt,
398 TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
401 triangle_mesh_parameters.internal_closed_curve_pt() = hole_pt;
404 double uniform_element_area=0.2;
405 triangle_mesh_parameters.element_area() = uniform_element_area;
409 TriangleMesh<ELEMENT>(triangle_mesh_parameters);
412 Problem::mesh_pt()=My_mesh_pt;
415 complete_problem_setup();
419 sprintf(filename,
"RESLT/trace.dat");
420 Trace_file.open(filename);
423 oomph_info <<
"Number of equations: "
424 << this->assign_eqn_numbers() << std::endl;
435 template<
class ELEMENT>
442 unsigned nbound=My_mesh_pt->nboundary();
443 for(
unsigned ibound=0;ibound<nbound;ibound++)
445 unsigned num_nod=My_mesh_pt->nboundary_node(ibound);
446 for (
unsigned inod=0;inod<num_nod;inod++)
449 Node* nod_pt=My_mesh_pt->boundary_node_pt(ibound,inod);
458 unsigned n_element = My_mesh_pt->nelement();
459 for(
unsigned e=0;e<n_element;e++)
462 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(My_mesh_pt->element_pt(e));
470 apply_boundary_conditions();
479 template<
class ELEMENT>
484 unsigned nbound=this->My_mesh_pt->nboundary();
485 for(
unsigned ibound=0;ibound<nbound;ibound++)
487 unsigned num_nod=this->My_mesh_pt->nboundary_node(ibound);
488 for (
unsigned inod=0;inod<num_nod;inod++)
491 Node* nod_pt=this->My_mesh_pt->boundary_node_pt(ibound,inod);
503 nod_pt->set_value(0,u[0]);
513 template<
class ELEMENT>
515 std::string& comment)
524 sprintf(filename,
"RESLT/soln%i.dat",Doc_info.number());
525 some_file.open(filename);
526 this->My_mesh_pt->output(some_file,npts);
527 some_file <<
"TEXT X = 22, Y = 92, CS=FRAME T = \""
528 << comment <<
"\"\n";
533 sprintf(filename,
"RESLT/exact_soln%i.dat",Doc_info.number());
534 some_file.open(filename);
540 sprintf(filename,
"RESLT/boundaries%i.dat",Doc_info.number());
541 some_file.open(filename);
542 My_mesh_pt->output_boundaries(some_file);
548 double error,norm,dummy_error,zero_norm;
549 sprintf(filename,
"RESLT/error%i.dat",Doc_info.number());
550 some_file.open(filename);
555 dummy_error,zero_norm);
559 oomph_info <<
"\nNorm of error : " << sqrt(error) << std::endl;
560 oomph_info <<
"Norm of exact solution: " << sqrt(norm) << std::endl;
561 oomph_info <<
"Norm of computed solution: " << sqrt(dummy_error) << std::endl;
562 Trace_file << sqrt(norm) <<
" " << sqrt(dummy_error) << std::endl;
573 int main(
int argc,
char **argv)
576 CommandLineArgs::setup(argc,argv);
582 CommandLineArgs::specify_command_line_flag(
"--validation");
585 CommandLineArgs::parse_and_assign();
588 CommandLineArgs::doc_specified_flags();
596 if (CommandLineArgs::command_line_flag_has_been_set(
"--validation"))
600 for (
unsigned i=0;i<nstep;i++)
604 problem.newton_solve();
608 std::stringstream comment_stream;
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
TriangleMesh< ELEMENT > * My_mesh_pt
Pointers to specific mesh.
UnstructuredPoissonProblem()
Constructor.
void actions_before_newton_solve()
Update the problem specs before solve: Re-apply boundary conditons.
void complete_problem_setup()
Helper function to (re-)set boundary condition and complete the build of all elements.
void actions_after_newton_solve()
Update after solve (empty)
void doc_solution(const std::string &comment="")
Doc the solution.
void apply_boundary_conditions()
Helper function to apply boundary conditions.
~UnstructuredPoissonProblem()
Destructor.
int main(int argc, char **argv)
Driver code for demo of inline triangle mesh generation.
Namespace for exact solution for Poisson equation with "sharp step".
void zero(const Vector< double > &x, Vector< double > &u)
Zero function – used to compute norm of the computed solution by computing the norm of the error when...
double TanPhi
Parameter for angle Phi of "step".
void get_source(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.