line_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 // oomph-lib headers
31 #include "map_matrix.h"
32 #include "line_mesh.h"
33 
34 namespace oomph
35 {
36 
37  //=======================================================================
38  /// Set up lookup schemes which establish which elements are located
39  /// next to which boundaries (doc to outfile if it's open)
40  //=======================================================================
41  void LineMeshBase::setup_boundary_element_info(std::ostream &outfile)
42  {
43  // Initialise documentation flag
44  bool doc = false;
45 
46  // Set this to true if an open file has been passed to the function
47  if(outfile) { doc = true; }
48 
49  // Determine number of boundaries in mesh
50  const unsigned n_bound = nboundary();
51 
52  // Wipe/allocate storage for arrays
53  Boundary_element_pt.clear();
54  Face_index_at_boundary.clear();
55  Boundary_element_pt.resize(n_bound);
56  Face_index_at_boundary.resize(n_bound);
57 
58  // Matrix map for working out the fixed local coord for elements on boundary
60  boundary_identifier;
61 
62  // Determine number of elements in the mesh
63  const unsigned n_element = nelement();
64 
65  // Loop over elements
66  for(unsigned e=0;e<n_element;e++)
67  {
68  // Get pointer to element
70 
71  // Output information to output file
72  if(doc) { outfile << "Element: " << e << " " << fe_pt << std::endl; }
73 
74  // Only include 1D elements! Some meshes contain interface elements too.
75  if(fe_pt->dim()==1)
76  {
77  // Loop over the element's nodes and find out which boundaries they're on
78  // ----------------------------------------------------------------------
79 
80  // Determine number of nodes in the element
81  const unsigned n_node = fe_pt->nnode_1d();
82 
83  // Loop over nodes in order
84  for (unsigned n=0;n<n_node;n++)
85  {
86  // Allocate storage for pointer to set of boundaries that node lives on
87  std::set<unsigned>* boundaries_pt = 0;
88 
89  // Get pointer to vector of boundaries that this node lives on
90  fe_pt->node_pt(n)->get_boundaries_pt(boundaries_pt);
91 
92  // If the node lives on some boundaries....
93  if (boundaries_pt!=0)
94  {
95  // Determine number of boundaries which node lives on
96  const unsigned n_bound_node_lives_on = (*boundaries_pt).size();
97 
98  // Throw an error if the node lives on more than one boundary
99  if(n_bound_node_lives_on>1)
100  {
101  std::string error_message =
102  "In a 1D mesh a node shouldn't be able to live on more than\n";
103  error_message +=
104  "one boundary, yet this node claims to.";
105 
106  throw OomphLibError(error_message,
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110  // If the node lives on just one boundary
111  else if(n_bound_node_lives_on==1)
112  {
113  // Determine which boundary the node lives on
114  const std::set<unsigned>::iterator boundary
115  = boundaries_pt->begin();
116 
117  // In 1D if an element has any nodes on a boundary then it must
118  // be a boundary element. This means that (unlike in the 2D and
119  // 3D cases) we can immediately add this element to permanent
120  // storage.
121  Boundary_element_pt[*boundary].push_back(fe_pt);
122 
123  // Record information required for FaceElements.
124  // `Face_index_at_boundary' = -/+1 for nodes on the left/right
125  // boundary. This allows us to decide which edge of the element
126  // coincides with the boundary since the line element must have
127  // precisely one vertex node on the boundary.
128 
129  // Are we at the left-hand vertex node? (left face)
130  if(n==0)
131  {
132  Face_index_at_boundary[*boundary].push_back(-1);
133  }
134 
135  // Are we at the right-hand vertex node? (right face)
136  else if(n==n_node-1)
137  {
138  Face_index_at_boundary[*boundary].push_back(1);
139  }
140  }
141  } // End of if node lives on some boundaries
142 
143  } // End of loop over nodes
144  }
145  } // End of loop over elements
146 
147 #ifdef PARANOID
148 
149  // Check each boundary only has precisely one element next to it
150  // -------------------------------------------------------------
151  //Only if not distributed
152 #ifdef OOMPH_HAS_MPI
153  if(Comm_pt==0)
154 #endif
155  {
156  // Loop over boundaries
157  for(unsigned b=0;b<n_bound;b++)
158  {
159  // Evaluate number of elements adjacent to boundary b
160  const unsigned n_element = Boundary_element_pt[b].size();
161 
162  switch(n_element)
163  {
164  // Boundary b has no adjacent elements
165  case 0:
166  {
167  std::ostringstream error_stream;
168  error_stream << "Boundary " << b << " has no element adjacent to it\n";
169  throw OomphLibError(error_stream.str(),
170  OOMPH_CURRENT_FUNCTION,
171  OOMPH_EXCEPTION_LOCATION);
172  break;
173  }
174  // Boundary b has one adjacent element (this is good!)
175  case 1:
176  break;
177 
178  // Boundary b has more than one adjacent element
179  default:
180  {
181  std::ostringstream error_stream;
182  error_stream << "Boundary " << b << " has " << n_element
183  << " elements adjacent to it.\n"
184  << "This shouldn't occur in a 1D mesh.\n";
185  throw OomphLibError(error_stream.str(),
186  OOMPH_CURRENT_FUNCTION,
187  OOMPH_EXCEPTION_LOCATION);
188  break;
189  }
190  } // End of switch
191 
192  // Because each boundary should only have one element adjacent to it,
193  // each `Face_index_at_boundary[b]' should be of size one.
194 
195  const unsigned face_index_at_boundary_size
196  = Face_index_at_boundary[b].size();
197 
198  if(face_index_at_boundary_size != 1)
199  {
200  std::ostringstream error_stream;
201  error_stream
202  << "Face_index_at_boundary[" << b << "] has size"
203  << face_index_at_boundary_size
204  << " which does not make sense.\n"
205  << "In a 1D mesh its size should always be one since only\n"
206  << "one element can be adjacent to any particular boundary";
207  throw OomphLibError(error_stream.str(),
208  OOMPH_CURRENT_FUNCTION,
209  OOMPH_EXCEPTION_LOCATION);
210  }
211  } // End of loop over boundaries
212  }
213 #endif
214 
215  // Doc?
216  // ----
217  if(doc)
218  {
219  // Loop over boundaries
220  for(unsigned b=0;b<n_bound;b++)
221  {
222  const unsigned n_element = Boundary_element_pt[b].size();
223  outfile << "Boundary: " << b
224  << " is adjacent to " << n_element
225  << " elements" << std::endl;
226 
227  // Loop over elements on given boundary
228  for(unsigned e=0;e<n_element;e++)
229  {
230  FiniteElement* fe_pt = Boundary_element_pt[b][e];
231  outfile << "Boundary element:" << fe_pt
232  << " Face index on boundary is "
233  << Face_index_at_boundary[b][e] << std::endl;
234  }
235  }
236  } // End of if(doc)
237 
238 
239  // Lookup scheme has now been set up
241 
242  } // End of setup_boundary_element_info()
243 
244 } // End of namespace
void setup_boundary_element_info()
Definition: line_mesh.h:83
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
A general Finite Element class.
Definition: elements.h:1271
e
Definition: cfortran.h:575
bool Lookup_for_elements_next_boundary_is_setup
Definition: mesh.h:98
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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
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.
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
Definition: mesh.h:130