mesh_from_inline_triangle.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 #include <fenv.h>
27 
28 //Generic routines
29 #include "generic.h"
30 
31 // The equations
32 #include "poisson.h"
33 
34 // The mesh
35 #include "meshes/triangle_mesh.h"
36 
37 using namespace std;
38 using namespace oomph;
39 
40 
41 
42 //===== start_of_namespace=============================================
43 /// Namespace for exact solution for Poisson equation with "sharp step"
44 //=====================================================================
46 {
47 
48  /// Parameter for steepness of "step"
49  double Alpha=5.0;
50 
51  /// Parameter for angle Phi of "step"
52  double TanPhi=0.0;
53 
54  /// Exact solution as a Vector
55  void get_exact_u(const Vector<double>& x, Vector<double>& u)
56  {
57  u[0]=tanh(1.0-Alpha*(TanPhi*x[0]-x[1]));
58  }
59 
60  /// Source function required to make the solution above an exact solution
61  void get_source(const Vector<double>& x, double& source)
62  {
63  source = 2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
64  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*
65  Alpha*Alpha*TanPhi*TanPhi+2.0*tanh(-1.0+Alpha*(TanPhi*x[0]-x[1]))*
66  (1.0-pow(tanh(-1.0+Alpha*(TanPhi*x[0]-x[1])),2.0))*Alpha*Alpha;
67  }
68 
69 
70  /// Zero function -- used to compute norm of the computed solution by
71  /// computing the norm of the error when compared against this.
72  void zero(const Vector<double>& x, Vector<double>& u)
73  {
74  u[0]=0.0;
75  }
76 
77 } // end of namespace
78 
79 
80 
81 /// ////////////////////////////////////////////////////////
82 /// ////////////////////////////////////////////////////////
83 /// ////////////////////////////////////////////////////////
84 
85 
86 //==start_of_problem_class============================================
87 /// Class definition
88 //====================================================================
89 template<class ELEMENT>
90 class UnstructuredPoissonProblem : public virtual Problem
91 {
92 
93 public:
94 
95  /// Constructor
97 
98  /// Destructor
100 
101  /// Actions before adapt. Empty
103 
104  /// Actions after adapt:
105  /// Setup the problem again -- remember that the mesh has been
106  /// completely rebuilt and its element's don't have any
107  /// pointers to source fcts etc. yet
109  {
110  complete_problem_setup();
111  }
112 
113  /// Update after solve (empty)
115 
116  /// Update the problem specs before solve: Re-apply boundary conditons
118  {
119  apply_boundary_conditions();
120  }
121 
122  /// Doc the solution
123  void doc_solution(const std::string& comment="");
124 
125 
126 private:
127 
128  /// Doc info object for labeling output
129  DocInfo Doc_info;
130 
131  /// Helper function to apply boundary conditions
132  void apply_boundary_conditions();
133 
134  /// Helper function to (re-)set boundary condition
135  /// and complete the build of all elements
136  void complete_problem_setup();
137 
138  /// Pointers to specific mesh
139  RefineableTriangleMesh<ELEMENT>* My_mesh_pt;
140 
141  /// Trace file to document norm of solution
142  ofstream Trace_file;
143 
144 }; // end_of_problem_class
145 
146 
147 
148 
149 
150 //==start_constructor=====================================================
151 /// Constructor
152 //========================================================================
153 template<class ELEMENT>
155 {
156  // Intrinsic coordinate along GeomObject
157  Vector<double> zeta(1);
158 
159  // Position vector on GeomObject
160  Vector<double> posn(2);
161 
162  // Ellipse defining the outer boundary
163  double x_center = 0.0;
164  double y_center = 0.0;
165  double A = 1.0;
166  double B = 1.0;
167  Ellipse * outer_boundary_ellipse_pt = new Ellipse(A,B);
168 
169  // Pointer to the closed curve that defines the outer boundary
170  TriangleMeshClosedCurve* closed_curve_pt=0;
171 
172  // Build outer boundary as Polygon?
173  //---------------------------------
174  bool polygon_for_outer_boundary=false;
175 #ifdef OUTER_POLYGON
176  polygon_for_outer_boundary=true;
177 #endif
178  if (polygon_for_outer_boundary)
179  {
180  // Number of segments that make up the boundary
181  unsigned n_seg = 5;
182  double unit_zeta = 0.5*MathematicalConstants::Pi/double(n_seg);
183 
184  // The boundary is bounded by two distinct boundaries, each
185  // represented by its own polyline
186  Vector<TriangleMeshCurveSection*> boundary_polyline_pt(2);
187 
188  // Vertex coordinates on boundary
189  Vector<Vector<double> > bound_coords(n_seg+1);
190 
191  // First part of the boundary
192  //---------------------------
193  for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
194  {
195  // Resize the vector
196  bound_coords[ipoint].resize(2);
197 
198  // Get the coordinates
199  zeta[0]=unit_zeta*double(ipoint);
200  outer_boundary_ellipse_pt->position(zeta,posn);
201  bound_coords[ipoint][0]=posn[0]+x_center;
202  bound_coords[ipoint][1]=posn[1]+y_center;
203  }
204 
205  // Build the 1st boundary polyline
206  unsigned boundary_id=0;
207  boundary_polyline_pt[0]=new TriangleMeshPolyLine(bound_coords,boundary_id);
208 
209  // Second part of the boundary
210  //----------------------------
211  unit_zeta*=3.0;
212  for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
213  {
214  // Resize the vector
215  bound_coords[ipoint].resize(2);
216 
217  // Get the coordinates
218  zeta[0]=(unit_zeta*double(ipoint))+0.5*MathematicalConstants::Pi;
219  outer_boundary_ellipse_pt->position(zeta,posn);
220  bound_coords[ipoint][0]=posn[0]+x_center;
221  bound_coords[ipoint][1]=posn[1]+y_center;
222  }
223 
224  // Build the 2nd boundary polyline
225  boundary_id=1;
226  boundary_polyline_pt[1]=new TriangleMeshPolyLine(bound_coords,boundary_id);
227 
228 
229  // Create the triangle mesh polygon for outer boundary
230  //----------------------------------------------------
231  TriangleMeshPolygon *outer_polygon =
232  new TriangleMeshPolygon(boundary_polyline_pt);
233 
234  // Enable redistribution of polylines
235  outer_polygon->
236  enable_redistribution_of_segments_between_polylines();
237 
238  // Set the pointer
239  closed_curve_pt = outer_polygon;
240 
241  }
242  // Build outer boundary as curvilinear
243  //------------------------------------
244  else
245  {
246 
247  // Provide storage for pointers to the two parts of the curvilinear boundary
248  Vector<TriangleMeshCurveSection*> outer_curvilinear_boundary_pt(2);
249 
250  // First bit
251  //----------
252  double zeta_start=0.0;
253  double zeta_end=MathematicalConstants::Pi;
254  unsigned nsegment=5;
255  unsigned boundary_id=0;
256  outer_curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
257  outer_boundary_ellipse_pt,zeta_start,zeta_end,nsegment,boundary_id);
258 
259  // Second bit
260  //-----------
261  zeta_start=MathematicalConstants::Pi;
262  zeta_end=2.0*MathematicalConstants::Pi;
263  nsegment=8;
264  boundary_id=1;
265  outer_curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
266  outer_boundary_ellipse_pt,zeta_start,zeta_end,nsegment,boundary_id);
267 
268  // Combine to curvilinear boundary and define the
269  //--------------------------------
270  // outer boundary
271  //--------------------------------
272  closed_curve_pt=
273  new TriangleMeshClosedCurve(outer_curvilinear_boundary_pt);
274 
275  }
276 
277 
278  // Now build the holes
279  //====================
280  Vector<TriangleMeshClosedCurve*> hole_pt(2);
281 
282  // Build polygonal hole
283  //=====================
284 
285  // Build first hole: A circle
286  x_center = 0.0;
287  y_center = 0.5;
288  A = 0.1;
289  B = 0.1;
290  Ellipse* polygon_ellipse_pt=new Ellipse(A,B);
291 
292  // Number of segments defining upper and lower half of the hole
293  unsigned n_seg = 6;
294  double unit_zeta = MathematicalConstants::Pi/double(n_seg);
295 
296  // This hole is bounded by two distinct boundaries, each
297  // represented by its own polyline
298  Vector<TriangleMeshCurveSection*> hole_polyline_pt(2);
299 
300 
301  // First boundary of polygonal hole
302  //---------------------------------
303 
304  // Vertex coordinates
305  Vector<Vector<double> > bound_hole(n_seg+1);
306  for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
307  {
308  // Resize the vector
309  bound_hole[ipoint].resize(2);
310 
311  // Get the coordinates
312  zeta[0]=unit_zeta*double(ipoint);
313  polygon_ellipse_pt->position(zeta,posn);
314  bound_hole[ipoint][0]=posn[0]+x_center;
315  bound_hole[ipoint][1]=posn[1]+y_center;
316  }
317 
318  // Specify the hole boundary id
319  unsigned boundary_id=2;
320 
321  // Build the 1st hole polyline
322  hole_polyline_pt[0] = new TriangleMeshPolyLine(bound_hole,boundary_id);
323 
324 
325  // Second boundary of polygonal hole
326  //----------------------------------
327  for(unsigned ipoint=0; ipoint<n_seg+1;ipoint++)
328  {
329  // Resize the vector
330  bound_hole[ipoint].resize(2);
331 
332  // Get the coordinates
333  zeta[0]=(unit_zeta*double(ipoint))+MathematicalConstants::Pi;
334  polygon_ellipse_pt->position(zeta,posn);
335  bound_hole[ipoint][0]=posn[0]+x_center;
336  bound_hole[ipoint][1]=posn[1]+y_center;
337  }
338 
339  // Specify the hole boundary id
340  boundary_id=3;
341 
342  // Build the 2nd hole polyline
343  hole_polyline_pt[1] = new TriangleMeshPolyLine(bound_hole,boundary_id);
344 
345 
346  // Build the polygonal hole
347  //-------------------------
348 
349  // Inner hole center coordinates
350  Vector<double> hole_center(2);
351  hole_center[0]=x_center;
352  hole_center[1]=y_center;
353 
354  hole_pt[0] = new TriangleMeshPolygon(hole_polyline_pt, hole_center);
355 
356 
357  // Build curvilinear hole
358  //======================
359 
360  // Build second hole: Another ellipse
361  A = 0.2;
362  B = 0.1;
363  Ellipse* ellipse_pt=new Ellipse(A,B);
364 
365  // Build the two parts of the curvilinear boundary
366  Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt(2);
367 
368 
369  // First part of curvilinear boundary
370  //-----------------------------------
371  double zeta_start=0.0;
372  double zeta_end=MathematicalConstants::Pi;
373  unsigned nsegment=10;
374  boundary_id=4;
375  curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
376  ellipse_pt,zeta_start,zeta_end,
377  nsegment,boundary_id);
378 
379  // Second part of curvilinear boundary
380  //-------------------------------------
381  zeta_start=MathematicalConstants::Pi;
382  zeta_end=2.0*MathematicalConstants::Pi;
383  nsegment=15;
384  boundary_id=5;
385  curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
386  ellipse_pt,zeta_start,zeta_end,
387  nsegment,boundary_id);
388 
389 
390  // Combine to hole
391  //----------------
392  Vector<double> hole_coords(2);
393  hole_coords[0]=0.0;
394  hole_coords[1]=0.0;
395  Vector<TriangleMeshClosedCurve*> curvilinear_hole_pt(1);
396  hole_pt[1]=
397  new TriangleMeshClosedCurve(curvilinear_boundary_pt,
398  hole_coords);
399 
400  // Uncomment this as an exercise to observe how a
401  // layer of fine elements get left behind near the boundary
402  // once the tanh step has swept past:
403 
404  // closed_curve_pt->disable_polyline_refinement();
405  // closed_curve_pt->disable_polyline_unrefinement();
406 
407  // Now build the mesh
408  //===================
409 
410  // Use the TriangleMeshParameters object for helping on the manage of the
411  // TriangleMesh parameters
412  TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
413 
414  // Specify the closed curve using the TriangleMeshParameters object
415  triangle_mesh_parameters.internal_closed_curve_pt() = hole_pt;
416 
417  // Specify the maximum area element
418  double uniform_element_area=0.2;
419  triangle_mesh_parameters.element_area() = uniform_element_area;
420 
421  // Create the mesh
422  My_mesh_pt=new
423  RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters);
424 
425  // Store as the problem's one and only mesh
426  Problem::mesh_pt()=My_mesh_pt;
427 
428  // Set error estimator for bulk mesh
429  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
430  My_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
431 
432  // Set element size limits
433  My_mesh_pt->max_element_size()=0.2;
434  My_mesh_pt->min_element_size()=0.002;
435 
436  // Set boundary condition and complete the build of all elements
437  complete_problem_setup();
438 
439  // Open trace file
440  char filename[100];
441  sprintf(filename,"RESLT/trace.dat");
442  Trace_file.open(filename);
443 
444  // Setup equation numbering scheme
445  oomph_info <<"Number of equations: "
446  << this->assign_eqn_numbers() << std::endl;
447 
448 } // end_of_constructor
449 
450 
451 
452 
453 //==start_of_complete======================================================
454  /// Set boundary condition exactly, and complete the build of
455  /// all elements
456 //========================================================================
457 template<class ELEMENT>
459 {
460 
461  // Set the boundary conditions for problem: All nodes are
462  // free by default -- just pin the ones that have Dirichlet conditions
463  // here.
464  unsigned nbound=My_mesh_pt->nboundary();
465  for(unsigned ibound=0;ibound<nbound;ibound++)
466  {
467  unsigned num_nod=My_mesh_pt->nboundary_node(ibound);
468  for (unsigned inod=0;inod<num_nod;inod++)
469  {
470  // Get node
471  Node* nod_pt=My_mesh_pt->boundary_node_pt(ibound,inod);
472 
473  // Pin one-and-only unknown value
474  nod_pt->pin(0);
475  }
476  } // end loop over boundaries
477 
478 
479  // Complete the build of all elements so they are fully functional
480  unsigned n_element = My_mesh_pt->nelement();
481  for(unsigned e=0;e<n_element;e++)
482  {
483  // Upcast from GeneralisedElement to the present element
484  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(My_mesh_pt->element_pt(e));
485 
486  //Set the source function pointer
487  el_pt->source_fct_pt() = &TanhSolnForPoisson::get_source;
488  }
489 
490  // Re-apply Dirichlet boundary conditions (projection ignores
491  // boundary conditions!)
492  apply_boundary_conditions();
493 }
494 
495 
496 
497 
498 //==start_of_apply_bc=====================================================
499  /// Helper function to apply boundary conditions
500 //========================================================================
501 template<class ELEMENT>
503 {
504 
505  // Loop over all boundary nodes
506  unsigned nbound=this->My_mesh_pt->nboundary();
507  for(unsigned ibound=0;ibound<nbound;ibound++)
508  {
509  unsigned num_nod=this->My_mesh_pt->nboundary_node(ibound);
510  for (unsigned inod=0;inod<num_nod;inod++)
511  {
512  // Get node
513  Node* nod_pt=this->My_mesh_pt->boundary_node_pt(ibound,inod);
514 
515  // Extract nodal coordinates from node:
516  Vector<double> x(2);
517  x[0]=nod_pt->x(0);
518  x[1]=nod_pt->x(1);
519 
520  // Compute the value of the exact solution at the nodal point
521  Vector<double> u(1);
523 
524  // Assign the value to the one (and only) nodal value at this node
525  nod_pt->set_value(0,u[0]);
526  }
527  }
528 
529 } // end set bc
530 
531 
532 //==start_of_doc_solution=================================================
533 /// Doc the solution
534 //========================================================================
535 template<class ELEMENT>
537  std::string& comment)
538 {
539  ofstream some_file;
540  char filename[100];
541 
542  // Number of plot points
543  unsigned npts;
544  npts=5;
545 
546  sprintf(filename,"RESLT/soln%i.dat",Doc_info.number());
547  some_file.open(filename);
548  this->My_mesh_pt->output(some_file,npts);
549  some_file << "TEXT X = 22, Y = 92, CS=FRAME T = \""
550  << comment << "\"\n";
551  some_file.close();
552 
553  // Output exact solution
554  //----------------------
555  sprintf(filename,"RESLT/exact_soln%i.dat",Doc_info.number());
556  some_file.open(filename);
557  My_mesh_pt->output_fct(some_file,npts,TanhSolnForPoisson::get_exact_u);
558  some_file.close();
559 
560  // Output boundaries
561  //------------------
562  sprintf(filename,"RESLT/boundaries%i.dat",Doc_info.number());
563  some_file.open(filename);
564  My_mesh_pt->output_boundaries(some_file);
565  some_file.close();
566 
567 
568  // Doc error and return of the square of the L2 error
569  //---------------------------------------------------
570  double error,norm,dummy_error,zero_norm;
571  sprintf(filename,"RESLT/error%i.dat",Doc_info.number());
572  some_file.open(filename);
573  My_mesh_pt->compute_error(some_file,TanhSolnForPoisson::get_exact_u,
574  error,norm);
575 
576  My_mesh_pt->compute_error(some_file,TanhSolnForPoisson::zero,
577  dummy_error,zero_norm);
578  some_file.close();
579 
580  // Doc L2 error and norm of solution
581  oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;
582  oomph_info << "Norm of exact solution: " << sqrt(norm) << std::endl;
583  oomph_info << "Norm of computed solution: " << sqrt(dummy_error) << std::endl;
584  Trace_file << sqrt(norm) << " " << sqrt(dummy_error) << std::endl;
585 
586  // Increment the doc_info number
587  Doc_info.number()++;
588 
589 } // end of doc
590 
591 
592 //=======start_of_main========================================
593 /// Driver code for demo of inline triangle mesh generation
594 //============================================================
595 int main(int argc, char **argv)
596 {
597  // Store command line arguments
598  CommandLineArgs::setup(argc,argv);
599 
600  // Define possible command line arguments and parse the ones that
601  // were actually specified
602 
603  // Validation?
604  CommandLineArgs::specify_command_line_flag("--validation");
605 
606  // Parse command line
607  CommandLineArgs::parse_and_assign();
608 
609  // Doc what has actually been specified on the command line
610  CommandLineArgs::doc_specified_flags();
611 
612  // Create problem
614  problem;
615 
616 // // Solve, adapt and doc manually
617 // unsigned nadapt=4;
618 // for (unsigned i=0;i<nadapt;i++)
619 // {
620 // problem.newton_solve();
621 // std::stringstream comment_stream;
622 // comment_stream << "Solution after " << i << " manual adaptations";
623 // problem.doc_solution(comment_stream.str());
624 // if (i!=(nadapt-1)) problem.adapt();
625 // }
626 
627 
628  // Doc the initial mesh
629  //=====================
630  std::stringstream comment_stream;
631  comment_stream << "Initial mesh ";
632  problem.doc_solution(comment_stream.str());
633 
634 
635  // Loop over orientation of step
636  //==============================
637  unsigned nstep=5;
638  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
639  {
640  nstep=2;
641  }
642  for (unsigned i=0;i<nstep;i++)
643  {
644  // Solve with spatial adaptation
645  //==============================
646  unsigned max_adapt=3;
647  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
648  {
649  max_adapt=1;
650  }
651  problem.newton_solve(max_adapt);
652 
653  // Doc the solution
654  //=================
655  std::stringstream comment_stream;
656  comment_stream << "Solution for tan(phi) = " << TanhSolnForPoisson::TanPhi;
657  problem.doc_solution(comment_stream.str());
658 
659  // Rotate orientation of solution
661  }
662 
663 } //End of main
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
void actions_before_newton_solve()
Update the problem specs before solve: Re-apply boundary conditons.
DocInfo Doc_info
Doc info object for labeling output.
void complete_problem_setup()
Helper function to (re-)set boundary condition and complete the build of all elements.
void actions_after_adapt()
Actions after adapt: Setup the problem again – remember that the mesh has been completely rebuilt and...
void actions_after_newton_solve()
Update after solve (empty)
void doc_solution(const std::string &comment="")
Doc the solution.
ofstream Trace_file
Trace file to document norm of solution.
void actions_before_adapt()
Actions before adapt. Empty.
void apply_boundary_conditions()
Helper function to apply boundary conditions.
RefineableTriangleMesh< ELEMENT > * My_mesh_pt
Pointers to specific mesh.
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.