triangle_mesh.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// Version 1.0; svn revision $LastChangedRevision: 1216 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-08-06 15:23:30 +0100 (Sat, 06 Aug 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #include <algorithm>
32 #include "map_matrix.h"
34 #include "triangle_mesh.h"
35 
36 namespace oomph
37 {
38 #ifdef OOMPH_HAS_TRIANGLE_LIB
39 
40 //==============================================================
41 /// Dump the triangulateio structure to a dump file and
42 /// record boundary coordinates of boundary nodes
43 //==============================================================
44  void TriangleMeshBase::dump_triangulateio(std::ostream &dump_file)
45  {
47 
48 #ifdef OOMPH_HAS_MPI
49  // If the mesh is not distributed then process what follows
50  if (!this->is_mesh_distributed())
51  {
52 #endif // #ifdef OOMPH_HAS_MPI
53 
54  // Loop over all boundary nodes and dump out boundary coordinates
55  // if they exist
56  Vector<double> zeta(1);
57  unsigned nb=nboundary();
58  for (unsigned b=0;b<nb;b++)
59  {
61  {
62  dump_file << "1 # Boundary coordinate for boundary " << b
63  << " does exist\n";
64  unsigned nnod=nboundary_node(b);
65  dump_file << nnod << " # Number of dumped boundary nodes\n";
66  for (unsigned j=0;j<nnod;j++)
67  {
68  Node* nod_pt=boundary_node_pt(b,j);
69  nod_pt->get_coordinates_on_boundary(b,zeta);
70  dump_file << zeta[0] << std::endl;
71  }
72  dump_file << "-999 # Done boundary coords for boundary " << b << "\n";
73  }
74  else
75  {
76  dump_file << "0 # Boundary coordinate for boundary " << b
77  << " does not exist\n";
78 
79  }
80  }
81 
82 #ifdef OOMPH_HAS_MPI
83  }
84 #endif // #ifdef OOMPH_HAS_MPI
85 
86  }
87 
88 
89 //==============================================================
90 /// Regenerate the mesh from a dumped triangulateio file
91 /// and dumped boundary coordinates of boundary nodes
92 //==============================================================
94  &restart_file)
95  {
96 
97 
98 #ifdef PARANOID
99  // Record number of boundaries
100  unsigned nbound_old=nboundary();
101 #endif
102 
103  //Clear the existing triangulate io
105 
106  //Read the data into the file
108 
109  //Now remesh from the new data structure
111 
112 #ifdef OOMPH_HAS_MPI
113  // If the mesh is not distributed then process what follows
114  if (!this->is_mesh_distributed())
115  {
116 #endif // #ifdef OOMPH_HAS_MPI
117 
118 #ifdef PARANOID
119  // Record number of boundary nodes after remesh
120  unsigned nbound_new=nboundary();
121  if (nbound_new!=nbound_old)
122  {
123  std::ostringstream error_stream;
124  error_stream << "Number of boundaries before remesh from triangulateio, "
125  << nbound_new
126  << ",\ndoesn't match number boundaries afterwards, "
127  << nbound_old
128  << ". Have you messed \naround with boundary nodes in the "
129  << "derived mesh constructor (or after calling \nit)? If so,"
130  << " the dump/restart won't work as written at the moment.";
131  throw OomphLibError(error_stream.str(),
132  OOMPH_CURRENT_FUNCTION,
133  OOMPH_EXCEPTION_LOCATION);
134  }
135 #endif
136 
137 
138  // Loop over all boundary nodes and read boundary coordinates
139  // if they exist
140  Vector<double> zeta(1);
141  std::string input_string;
142  unsigned nb=nboundary();
143  for (unsigned b=0;b<nb;b++)
144  {
145  // Read line up to termination sign
146  getline(restart_file,input_string,'#');
147 
148  // Ignore rest of line
149  restart_file.ignore(80,'\n');
150 
151  // Did boundary coordinate exist?
152  const unsigned bound_coord_exists=atoi(input_string.c_str());
153  if (bound_coord_exists==1)
154  {
155  // Remember it!
157 
158  // Read line up to termination sign
159  getline(restart_file,input_string,'#');
160 
161  // Ignore rest of line
162  restart_file.ignore(80,'\n');
163 
164  // How many nodes did we dump?
165  const unsigned nnod_dumped=atoi(input_string.c_str());
166 
167  // Does it match?
168  unsigned nnod=nboundary_node(b);
169  if (nnod!=nnod_dumped)
170  {
171  std::ostringstream error_stream;
172  error_stream << "Number of dumped boundary nodes "
173  << nnod_dumped
174  << " doesn't match number of nodes on boundary "
175  << b << ": " << nnod << std::endl;
176  throw OomphLibError(error_stream.str(),
177  OOMPH_CURRENT_FUNCTION,
178  OOMPH_EXCEPTION_LOCATION);
179  }
180 
181  // Loop over all nodes
182  for (unsigned j=0;j<nnod;j++)
183  {
184  // Read line up to termination sign
185  getline(restart_file,input_string);
186 
187  // Boundary coordinate
188  zeta[0]=atof(input_string.c_str());
189 
190  // Set it
191  Node* nod_pt=boundary_node_pt(b,j);
192  nod_pt->set_coordinates_on_boundary(b,zeta);
193  }
194 
195  // Read line up to termination sign
196  getline(restart_file,input_string,'#');
197 
198  // Ignore rest of line
199  restart_file.ignore(80,'\n');
200 
201  // Have we reached the end?
202  const int check=atoi(input_string.c_str());
203  if (check!=-999)
204  {
205  std::ostringstream error_stream;
206  error_stream << "Haven't read all nodes on boundary "<< b
207  << std::endl;
208  throw OomphLibError(error_stream.str(),
209  OOMPH_CURRENT_FUNCTION,
210  OOMPH_EXCEPTION_LOCATION);
211  }
212  }
213  else
214  {
215  oomph_info << "Restart: Boundary coordinate for boundary " << b
216  << " does not exist.\n";
217 
218  }
219  }
220 #ifdef OOMPH_HAS_MPI
221  } // if (!this->is_mesh_distributed())
222 #endif // #ifdef OOMPH_HAS_MPI
223 
224  }
225 
226 
227 
228 
229 
230 //==============================================================
231 /// Write a Triangulateio_object file of the TriangulateIO object
232 /// String s is add to assign a different value for
233 /// input and/or output structure.
234 /// The function give the same result of the "report" function
235 /// included in the tricall.c, esternal_src.
236 //==============================================================
238  std::string &s)
239  {
240 
241  std::ofstream outfile;
242  char filename[100];
243 
244  sprintf(filename,"Triangulateio_object_%s.dat",s.c_str());
245  outfile.open(filename);
246  outfile <<"# Triangulateio object values:\n\n"<<std::endl;
247 
248  // Write points coordinates
249  if(triangle.numberofpoints!=0)
250  {
251  outfile <<"# Triangulateio number of points is:"
252  <<triangle.numberofpoints <<std::endl;
253  }
254  if(triangle.pointlist != NULL)
255  {
256  outfile <<"# Vertex coordinates are:"<<std::endl;
257  for(int k=0;k<triangle.numberofpoints*2;k+=2)
258  {
259  outfile << (k*0.5)+1 << " "
260  << triangle.pointlist[k] << " "
261  << triangle.pointlist[k+1] << std::endl;
262  }
263  }
264 
265  // Write points attribute list
266  if(triangle.numberofpointattributes!=0)
267  {
268  outfile <<"# Triangulateio number of points attributelist is:"
269  <<triangle.numberofpointattributes <<std::endl;
270  }
271  if(triangle.pointattributelist != NULL)
272  {
273 
274  outfile <<"# Vertex attribute are:"<<std::endl;
275  for(int k=0;k<triangle.numberofpointattributes;k++)
276  {
277  outfile << triangle.pointattributelist[k] << std::endl;
278  }
279  }
280 
281  // Write point markers list
282  if(triangle.pointmarkerlist != NULL)
283  {
284  outfile <<"# Vertex Markers are:"<<std::endl;
285  for(int k=0;k<triangle.numberofpoints;k++)
286  {
287  outfile << triangle.pointmarkerlist[k] << std::endl;
288  }
289  }
290 
291  // Write the 1.node file used by the showme function
292  std::ofstream nodefile;
293  char nodename[100];
294 
295  sprintf(nodename,"file_%s.1.node",s.c_str());
296  nodefile.open(nodename);
297  nodefile <<triangle.numberofpoints<<" 2 "
298  <<triangle.numberofpointattributes
299  <<" 0"<<std::endl;
300  for(int j=0;j<triangle.numberofpoints*2;j+=2)
301  {
302  nodefile <<(j/2)+1<<" " << triangle.pointlist[j] << " "
303  << triangle.pointlist[j+1] << std::endl;
304  }
305  nodefile.close();
306 
307 
308  // Write segments edge elements
309  if(triangle.numberofsegments!=0)
310  {
311  outfile <<"# The number of segments is:"
312  <<triangle.numberofsegments<<std::endl;
313  }
314  if(triangle.segmentlist != NULL)
315  {
316  outfile <<"# Segments are:"<<std::endl;
317  for(int k=0;k<triangle.numberofsegments*2;k+=2)
318  {
319  outfile << triangle.segmentlist[k] << " "
320  << triangle.segmentlist[k+1] <<std::endl;
321  }
322  }
323 
324  // Write segments markers list
325  if(triangle.segmentmarkerlist != NULL)
326  {
327  outfile <<"# Segments Markers are:"<<std::endl;
328  for(int k=0;k<triangle.numberofsegments;k++)
329  {
330  outfile << triangle.segmentmarkerlist[k] << std::endl;
331  }
332  }
333 
334  // Write regions
335  if(triangle.numberofregions!=0)
336  {
337  outfile <<"# The number of region is:"
338  <<triangle.numberofregions<<std::endl;
339  }
340 
341  // Write holes
342  if(triangle.numberofholes!=0)
343  {
344  outfile <<"# The number of holes is:"
345  <<triangle.numberofholes<<std::endl;
346  }
347  if(triangle.holelist != NULL)
348  {
349  outfile <<"# Holes are:"<<std::endl;
350  for(int k=0;k<triangle.numberofholes*2;k+=2)
351  {
352  outfile << triangle.holelist[k] << " "
353  << triangle.holelist[k+1] <<std::endl;
354  }
355  }
356 
357  // Write triangles
358  if(triangle.numberoftriangles!=0)
359  {
360  outfile <<"# Triangulateio number of triangles:"
361  <<triangle.numberoftriangles <<std::endl;
362  }
363  if(triangle.numberofcorners!=0)
364  {
365  outfile <<"# Triangulateio number of corners:"
366  <<triangle.numberofcorners<<std::endl;
367  }
368  if(triangle.numberoftriangleattributes!=0)
369  {
370  outfile <<"# Triangulateio number of triangles attributes:"
371  <<triangle.numberoftriangleattributes <<std::endl;
372  }
373  if(triangle.trianglelist != NULL)
374  {
375 
376  outfile <<"# Traingles are:"<<std::endl;
377  for(int k=0;k<triangle.numberoftriangles*3;k+=3)
378  {
379  outfile << triangle.trianglelist[k] << " "
380  << triangle.trianglelist[k+1] << " "
381  << triangle.trianglelist[k+2] << std::endl;
382  }
383  }
384 
385  if(triangle.trianglearealist != NULL)
386  {
387  outfile <<"# Triangle's areas are:"<<std::endl;
388  for(int k=0;k<triangle.numberoftriangles;k++)
389  {
390  outfile << triangle.trianglearealist[k] << std::endl;
391  }
392  }
393 
394  if(triangle.trianglelist != NULL)
395  {
396 
397  // Write the 1.ele file used by the showme function
398  std::ofstream elefile;
399  char elename[100];
400 
401  sprintf(elename,"file_%s.1.ele",s.c_str());
402  elefile.open(elename);
403  elefile <<triangle.numberoftriangles<<" 3 0"<<std::endl;
404  for(int j=0;j<triangle.numberoftriangles*3;j+=3)
405  {
406  elefile <<(j/3)+1<<" "<< triangle.trianglelist[j] << " "
407  << triangle.trianglelist[j+1] << " "
408  << triangle.trianglelist[j+2] << std::endl;
409  }
410  elefile.close();
411  }
412 
413  outfile.close();
414  }
415 
416 #endif
417 
418 //================================================================
419 /// Setup lookup schemes which establish which elements are located
420 /// next to which boundaries (Doc to outfile if it's open).
421 //================================================================
423  &outfile)
424  {
425  // Should we document the output here
426  bool doc = false;
427 
428  if(outfile) doc = true;
429 
430  // Number of boundaries
431  unsigned nbound=nboundary();
432 
433  // Wipe/allocate storage for arrays
434  Boundary_element_pt.clear();
435  Face_index_at_boundary.clear();
436  Boundary_element_pt.resize(nbound);
437  Face_index_at_boundary.resize(nbound);
438 
439  // Temporary vector of vectors of pointers to elements on the boundaries:
440  // This is a vector to ensure that order is strictly preserved
441  Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
442  vector_of_boundary_element_pt.resize(nbound);
443 
444  // Matrix map for working out the fixed face for elements on boundary
446  face_identifier;
447 
448  // Loop over elements
449  //-------------------
450  unsigned nel=nelement();
451 
452  // Get pointer to vector of boundaries that the
453  // node lives on
454  Vector<std::set<unsigned>*> boundaries_pt(3,0);
455 
456  // Data needed to deal with edges through the
457  // interior of the domain
458  std::map<Edge,unsigned > edge_count;
459  std::map<Edge,TriangleBoundaryHelper::BCInfo>
460  edge_bcinfo;
461  std::map<Edge,TriangleBoundaryHelper::BCInfo>
462  face_info;
464  Vector<unsigned> bonus(nbound);
465 
466  // When using internal boundaries, an edge can be related to more than
467  // one element (because of both sides of the internal boundaries)
468  std::map<Edge,Vector<TriangleBoundaryHelper::BCInfo> > edge_internal_bnd;
469 
470  for (unsigned e=0;e<nel;e++)
471  {
472  // Get pointer to element
474 
475  if (doc) {outfile << "Element: " << e << " " << fe_pt << std::endl;}
476 
477  // Only include 2D elements! Some meshes contain interface elements too.
478  if (fe_pt->dim()==2)
479  {
480  // Loop over the element's nodes and find out which boundaries they're on
481  // ----------------------------------------------------------------------
482 
483  //We need only loop over the corner nodes
484  for(unsigned i=0;i<3;i++)
485  {
486  fe_pt->node_pt(i)->get_boundaries_pt(boundaries_pt[i]);
487  }
488 
489  //Find the common boundaries of each edge
490  Vector<std::set<unsigned> > edge_boundary(3);
491 
492  //Edge 0 connects points 1 and 2
493  //-----------------------------
494 
495  if(boundaries_pt[1] && boundaries_pt[2])
496  {
497  // Create the corresponding edge
498  Edge edge0(fe_pt->node_pt(1),fe_pt->node_pt(2));
499 
500  // Update infos about this edge
502  info.Face_id=0;
503  info.FE_pt = fe_pt;
504 
505  std::set_intersection(boundaries_pt[1]->begin(),boundaries_pt[1]->end(),
506  boundaries_pt[2]->begin(),boundaries_pt[2]->end(),
507  std::insert_iterator<std::set<unsigned> >(
508  edge_boundary[0],edge_boundary[0].begin()));
509  std::set<unsigned>::iterator it0=edge_boundary[0].begin();
510 
511  // Edge does exist:
512  if ( edge_boundary[0].size() > 0 )
513  {
514  info.Boundary=*it0;
515 
516  // How many times this edge has been visited
517  edge_count[edge0]++;
518 
519  // Update edge_bcinfo
520  edge_bcinfo.insert(std::make_pair(edge0,info));
521 
522  // ... and also update the info associated with internal bnd
523  edge_internal_bnd[edge0].push_back(info);
524  }
525  }
526 
527  //Edge 1 connects points 0 and 2
528  //-----------------------------
529 
530  if(boundaries_pt[0] && boundaries_pt[2])
531  {
532  std::set_intersection(boundaries_pt[0]->begin(),boundaries_pt[0]->end(),
533  boundaries_pt[2]->begin(),boundaries_pt[2]->end(),
534  std::insert_iterator<std::set<unsigned> >(
535  edge_boundary[1],edge_boundary[1].begin()));
536 
537  // Create the corresponding edge
538  Edge edge1(fe_pt->node_pt(0),fe_pt->node_pt(2));
539 
540  // Update infos about this edge
542  info.Face_id=1;
543  info.FE_pt = fe_pt;
544  std::set<unsigned>::iterator it1=edge_boundary[1].begin();
545 
546  // Edge does exist:
547  if ( edge_boundary[1].size() > 0)
548  {
549  info.Boundary=*it1;
550 
551  // How many times this edge has been visited
552  edge_count[edge1]++;
553 
554  // Update edge_bcinfo
555  edge_bcinfo.insert(std::make_pair(edge1,info));
556 
557  // ... and also update the info associated with internal bnd
558  edge_internal_bnd[edge1].push_back(info);
559  }
560  }
561 
562  //Edge 2 connects points 0 and 1
563  //-----------------------------
564 
565  if(boundaries_pt[0] && boundaries_pt[1])
566  {
567  std::set_intersection(boundaries_pt[0]->begin(),boundaries_pt[0]->end(),
568  boundaries_pt[1]->begin(),boundaries_pt[1]->end(),
569  std::insert_iterator<std::set<unsigned> >(
570  edge_boundary[2],edge_boundary[2].begin()));
571 
572  // Create the corresponding edge
573  Edge edge2(fe_pt->node_pt(0),fe_pt->node_pt(1));
574 
575  // Update infos about this edge
577  info.Face_id=2;
578  info.FE_pt = fe_pt;
579  std::set<unsigned>::iterator it2=edge_boundary[2].begin();
580 
581  // Edge does exist:
582  if ( edge_boundary[2].size() > 0)
583  {
584  info.Boundary=*it2;
585 
586  // How many times this edge has been visited
587  edge_count[edge2]++;
588 
589  // Update edge_bcinfo
590  edge_bcinfo.insert(std::make_pair(edge2,info));
591 
592  // ... and also update the info associated with internal bnd
593  edge_internal_bnd[edge2].push_back(info);
594  }
595  }
596 
597 
598 #ifdef PARANOID
599 
600  // Check if edge is associated with multiple boundaries
601 
602  //We now know whether any edges lay on the boundaries
603  for(unsigned i=0;i<3;i++)
604  {
605  //How many boundaries are there
606  unsigned count = 0;
607 
608  //Loop over all the members of the set and add to the count
609  //and set the boundary
610  for(std::set<unsigned>::iterator it=edge_boundary[i].begin();
611  it!=edge_boundary[i].end();++it)
612  {
613  ++count;
614  }
615 
616  //If we're on more than one boundary, this is weird, so die
617  if(count > 1)
618  {
619  std::ostringstream error_stream;
620  error_stream << "Edge " << i << " is located on " <<
621  count << " boundaries.\n";
622  error_stream << "This is rather strange, so I'm going to die\n";
623  throw OomphLibError(
624  error_stream.str(),
625  OOMPH_CURRENT_FUNCTION,
626  OOMPH_EXCEPTION_LOCATION);
627  }
628  }
629 
630 #endif
631 
632  // Now we set the pointers to the boundary sets to zero
633  for(unsigned i=0;i<3;i++) {boundaries_pt[i] = 0;}
634 
635  }
636  } //end of loop over all elements
637 
638  // Loop over all edges that are located on a boundary
639  typedef std::map<Edge,TriangleBoundaryHelper::BCInfo>::iterator ITE;
640  for (ITE it=edge_bcinfo.begin();
641  it!=edge_bcinfo.end();
642  it++)
643  {
644  Edge current_edge = it->first;
645  unsigned bound=it->second.Boundary;
646 
647  // If the edge has been visited only once
648  if(edge_count[current_edge]==1)
649  {
650  // Count the edges that are on the same element and on the same boundary
651  face_count(static_cast<unsigned>(bound),it->second.FE_pt)=
652  face_count(static_cast<unsigned>(bound),it->second.FE_pt) + 1;
653 
654  //If such edges exist, let store the corresponding element
655  if( face_count(bound,it->second.FE_pt) > 1)
656  {
657  // Update edge's infos
659  info.Face_id=it->second.Face_id;
660  info.FE_pt = it->second.FE_pt;
661  info.Boundary=it->second.Boundary;
662 
663  // Add it to FIinfo, that stores infos of problematic elements
664  face_info.insert(std::make_pair(current_edge,info));
665 
666  //How many edges on which boundary have to be added
667  bonus[bound]++;
668  }
669  else
670  {
671  //Add element and face to the appropriate vectors
672  // Does the pointer already exits in the vector
674  std::find(vector_of_boundary_element_pt[
675  static_cast<unsigned>(bound)].begin(),
676  vector_of_boundary_element_pt[
677  static_cast<unsigned>(bound)].end(),
678  it->second.FE_pt);
679 
680  //Only insert if we have not found it (i.e. got to the end)
681  if(b_el_it == vector_of_boundary_element_pt[
682  static_cast<unsigned>(bound)].end())
683  {
684  vector_of_boundary_element_pt[static_cast<unsigned>(bound)].
685  push_back(it->second.FE_pt);
686  }
687 
688  //set_of_boundary_element_pt[static_cast<unsigned>(bound)].insert(
689  // it->second.FE_pt);
690  face_identifier(static_cast<unsigned>(bound),it->second.FE_pt) =
691  it->second.Face_id;
692  }
693 
694  }
695 
696  } //End of "adding-boundaries"-loop
697 
698 
699  // Now copy everything across into permanent arrays
700  //-------------------------------------------------
701 
702  // Loop over boundaries
703  for (unsigned i=0;i<nbound;i++)
704  {
705  // Number of elements on this boundary that have to be added
706  // in addition to other elements
707  unsigned bonus1=bonus[i];
708 
709  // Number of elements on this boundary
710  unsigned nel=vector_of_boundary_element_pt[i].size() + bonus1;
711 
712  // Allocate storage for the coordinate identifiers
713  Face_index_at_boundary[i].resize(nel);
714 
715  unsigned e_count=0;
717  for (IT it=vector_of_boundary_element_pt[i].begin();
718  it!=vector_of_boundary_element_pt[i].end();
719  it++)
720  {
721  // Recover pointer to element
722  FiniteElement* fe_pt=*it;
723 
724  // Add to permanent storage
725  Boundary_element_pt[i].push_back(fe_pt);
726 
727  Face_index_at_boundary[i][e_count] = face_identifier(i,fe_pt);
728 
729  // Increment counter
730  e_count++;
731 
732  }
733  // We add the elements that have two or more edges on this boundary
734  for (ITE itt= face_info.begin(); itt!= face_info.end(); itt++)
735  {
736  if (itt->second.Boundary==i)
737  {
738  // Add to permanent storage
739  Boundary_element_pt[i].push_back(itt->second.FE_pt);
740 
741  Face_index_at_boundary[i][e_count] = itt->second.Face_id;
742 
743  e_count++;
744  }
745 
746  }
747 
748  } //End of loop over boundaries
749 
750  // Doc?
751  //-----
752  if (doc)
753  {
754  // Loop over boundaries
755  for (unsigned i=0;i<nbound;i++)
756  {
757  unsigned nel=Boundary_element_pt[i].size();
758  outfile << "Boundary: " << i
759  << " is adjacent to " << nel
760  << " elements" << std::endl;
761 
762  // Loop over elements on given boundary
763  for (unsigned e=0;e<nel;e++)
764  {
766  outfile << "Boundary element:" << fe_pt
767  << " Face index of boundary is "
768  << Face_index_at_boundary[i][e] << std::endl;
769  }
770  }
771  }
772 
773  // Lookup scheme has now been setup yet
775  }
776 }
TriangulateIO Triangulateio
TriangulateIO representation of the mesh.
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.
double * pointattributelist
Pointer to list of point attributes.
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:497
int * pointmarkerlist
Pointer to list of point markers.
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:201
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1305
A general Finite Element class.
Definition: elements.h:1271
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
OomphInfo oomph_info
e
Definition: cfortran.h:575
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:98
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:809
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interface provides run-time error checking.
Definition: nodes.cc:2287
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1287
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:106
void clear_triangulateio(TriangulateIO &triangulate_io, const bool &clear_hole_data)
Clear TriangulateIO structure.
static char t char * s
Definition: cfortran.h:572
void remesh_from_triangulateio(std::istream &restart_file)
Regenerate the mesh from a dumped triangulateio file and dumped boundary coordinates of boundary node...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
FiniteElement * FE_pt
Pointer to bulk finite element.
Edge class.
Definition: mesh.h:2401
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,e)
Definition: mesh.h:102
void dump_triangulateio(TriangulateIO &triangle_io, std::ostream &dump_file)
Write all the triangulateio data to disk in a dump file that can then be used to restart simulations...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void write_triangulateio(TriangulateIO &triangulate_io, std::string &s)
Helper function. Write a TriangulateIO object file with all the triangulateio fields. String s is add to assign a different value for the input and/or output structure.
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
virtual void remesh_from_internal_triangulateio()
Virtual function that is used for specific remeshing from the triangulateio.
void read_triangulateio(std::istream &restart_file, TriangulateIO &triangle_io)
Read the triangulateio data from a dump file on disk, which can then be used to restart simulations...
void dump_triangulateio(std::ostream &dump_file)
Dump the triangulateio structure to a dump file and record boundary coordinates of boundary nodes...