quad_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: 1258 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-10-12 07:13:28 +0100 (Wed, 12 Oct 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 #include <algorithm>
31 #include "map_matrix.h"
32 #include "quad_mesh.h"
33 
34 
35 namespace oomph
36 {
37 
38 
39 //================================================================
40 /// Setup lookup schemes which establish which elements are located
41 /// next to which boundaries (Doc to outfile if it's open).
42 //================================================================
43 void QuadMeshBase::setup_boundary_element_info(std::ostream &outfile)
44 {
45 
46  bool doc=false;
47  if (outfile) doc=true;
48 
49  // Number of boundaries
50  unsigned nbound=nboundary();
51 
52  // Wipe/allocate storage for arrays
53  Boundary_element_pt.clear();
54  Face_index_at_boundary.clear();
55  Boundary_element_pt.resize(nbound);
56  Face_index_at_boundary.resize(nbound);
57 
58  // Temporary vector of vectors to pointers to elements on the boundaries:
59  // This is not a set to ensure UNIQUE ordering
60  Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
61  vector_of_boundary_element_pt.resize(nbound);
62 
63  // Matrix map for working out the fixed local coord for elements on boundary
65  boundary_identifier;
66 
67 
68  // Loop over elements
69  //-------------------
70  unsigned nel=nelement();
71  for (unsigned e=0;e<nel;e++)
72  {
73  // Get pointer to element
75 
76  if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
77 
78  // Only include 2D elements! Some meshes contain interface elements too.
79  if (fe_pt->dim()==2)
80  {
81  // Loop over the element's nodes and find out which boundaries they're on
82  // ----------------------------------------------------------------------
83  unsigned nnode_1d=fe_pt->nnode_1d();
84 
85  // Loop over nodes in order
86  for (unsigned i0=0;i0<nnode_1d;i0++)
87  {
88  for (unsigned i1=0;i1<nnode_1d;i1++)
89  {
90  // Local node number
91  unsigned j=i0+i1*nnode_1d;
92 
93  // Get pointer to vector of boundaries that this
94  // node lives on
95  std::set<unsigned>* boundaries_pt=0;
96  fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
97 
98  // If the node lives on some boundaries....
99  if (boundaries_pt!=0)
100  {
101  // Loop over boundaries
102  //unsigned nbound=(*boundaries_pt).size();
103  for (std::set<unsigned>::iterator it=boundaries_pt->begin();
104  it != boundaries_pt->end();++it)
105  {
106  // Add pointer to finite element to vector for the appropriate
107  // boundary
108 
109  // Does the pointer already exits in the vector
111  std::find(vector_of_boundary_element_pt[*it].begin(),
112  vector_of_boundary_element_pt[*it].end(),
113  fe_pt);
114  //Only insert if we have not found it (i.e. got to the end)
115  if(b_el_it == vector_of_boundary_element_pt[*it].end())
116  {
117  vector_of_boundary_element_pt[*it].push_back(fe_pt);
118  }
119 
120  // For the current element/boundary combination, create
121  // a vector that stores an indicator which element boundaries
122  // the node is located (boundary_identifier=-/+1 for nodes
123  // on the left/right boundary; boundary_identifier=-/+2 for nodes
124  // on the lower/upper boundary. We determine these indices
125  // for all corner nodes of the element and add them to a vector
126  // to a vector. This allows us to decide which face of the element
127  // coincides with the boundary since the (quad!) element must
128  // have exactly two corner nodes on the boundary.
129  if (boundary_identifier(*it,fe_pt)==0)
130  {
131  boundary_identifier(*it,fe_pt)=new Vector<int>;
132  }
133 
134  // Are we at a corner node?
135  if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1)))
136  {
137  // Create index to represent position relative to s_0
138  (*boundary_identifier(*it,fe_pt)).
139  push_back(1*(2*i0/(nnode_1d-1)-1));
140 
141  // Create index to represent position relative to s_1
142  (*boundary_identifier(*it,fe_pt)).
143  push_back(2*(2*i1/(nnode_1d-1)-1));
144  }
145 
146  }
147  }
148  //else
149  // {
150  // oomph_info << "...does not live on any boundaries " << std::endl;
151  // }
152 
153  }
154  }
155  }
156  //else
157  //{
158  //oomph_info << "Element " << e << " does not qualify" << std::endl;
159  //}
160  }
161 
162  // Now copy everything across into permanent arrays
163  //-------------------------------------------------
164 
165  // Note: vector_of_boundary_element_pt contains all elements
166  // that have (at least) one corner node on a boundary -- can't copy
167  // them across into Boundary_element_pt
168  // yet because some of them might have only one node on the
169  // the boundary in which case they don't qualify as
170  // boundary elements!
171 
172  // Loop over boundaries
173  //---------------------
174  for (unsigned i=0;i<nbound;i++)
175  {
176  // Number of elements on this boundary
177  // nel is unused, so I've commented it out - RWhite.
178 // unsigned nel=vector_of_boundary_element_pt[i].size();
179 
180  // Allocate storage for the face identifiers
181  //Face_index_at_boundary[i].resize(nel);
182 
183  // Loop over elements that have at least one corner node on this boundary
184  //-----------------------------------------------------------------------
185  //unsigned e_count=0;
187  for (IT it=vector_of_boundary_element_pt[i].begin();
188  it!=vector_of_boundary_element_pt[i].end();
189  it++)
190  {
191  // Recover pointer to element
192  FiniteElement* fe_pt=*it;
193 
194  // Initialise count for boundary identiers (-2,-1,1,2)
195  std::map<int,unsigned> count;
196 
197  // Loop over coordinates
198  for (unsigned ii=0;ii<2;ii++)
199  {
200  // Loop over upper/lower end of coordinates
201  for (int sign=-1;sign<3;sign+=2)
202  {
203  count[(ii+1)*sign]=0;
204  }
205  }
206 
207  // Loop over boundary indicators for this element/boundary
208  unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
209  for (unsigned j=0;j<n_indicators;j++)
210  {
211  count[(*boundary_identifier(i,fe_pt))[j] ]++;
212  }
213  delete boundary_identifier(i,fe_pt);
214 
215  // Determine the correct boundary indicator by checking that it
216  // occurs twice (since two corner nodes of the element's boundary
217  // need to be located on the domain boundary
218  int indicator=-10;
219 
220  //Check that we're finding exactly one boundary indicator
221  //bool found=false;
222 
223  // Loop over coordinates
224  for (unsigned ii=0;ii<2;ii++)
225  {
226  // Loop over upper/lower end of coordinates
227  for (int sign=-1;sign<3;sign+=2)
228  {
229  //If an index occurs twice then that face is on the boundary
230  //But we can have multiple faces on the same boundary id, so just add all the ones that we
231  //have
232  if (count[(ii+1)*sign]==2)
233  {
234  // Check that we haven't found multiple boundaries
235  /*if (found)
236  {
237  std::string error_message=
238  "Trouble: Multiple boundary identifiers!\n";
239  error_message +=
240  "Elements should only have at most 2 corner ";
241  error_message +=
242  "nodes on any one boundary.\n";
243 
244  throw OomphLibError(
245  error_message,
246  OOMPH_CURRENT_FUNCTION,
247  OOMPH_EXCEPTION_LOCATION);
248  }
249  found=true;*/
250  indicator=(ii+1)*sign;
251 
252  //Copy into the data structure
253  Boundary_element_pt[i].push_back(*it);
254  Face_index_at_boundary[i].push_back(indicator);
255  }
256  }
257  }
258 
259  // Element has exactly two corner nodes on boundary
260  /*if (found)
261  {
262  // Add to permanent storage
263  Boundary_element_pt[i].push_back(fe_pt);
264 
265  // Now convert boundary indicator into information required
266  // for FaceElements
267  switch (indicator)
268  {
269  //South face
270  case -2:
271 
272  // s_1 is fixed at -1.0:
273  Face_index_at_boundary[i][e_count] = -2;
274  break;
275 
276  //West face
277  case -1:
278 
279  // s_0 is fixed at -1.0:
280  Face_index_at_boundary[i][e_count] = -1;
281  break;
282 
283  //East face
284  case 1:
285 
286  // s_0 is fixed at 1.0:
287  Face_index_at_boundary[i][e_count] = 1;
288  break;
289 
290  //North face
291  case 2:
292 
293  // s_1 is fixed at 1.0:
294  Face_index_at_boundary[i][e_count] = 2;
295  break;
296 
297  default:
298 
299  throw OomphLibError("Never get here",
300  OOMPH_CURRENT_FUNCTION,
301  OOMPH_EXCEPTION_LOCATION);
302  }
303 
304  // Increment counter
305  e_count++;
306  }*/
307 
308  }
309  }
310 
311 
312 
313  // Doc?
314  //-----
315  if (doc)
316  {
317  // Loop over boundaries
318  for (unsigned i=0;i<nbound;i++)
319  {
320  unsigned nel=Boundary_element_pt[i].size();
321  outfile << "Boundary: " << i
322  << " is adjacent to " << nel
323  << " elements" << std::endl;
324 
325  // Loop over elements on given boundary
326  for (unsigned e=0;e<nel;e++)
327  {
329  outfile << "Boundary element:" << fe_pt
330  << " Face index on boundary is "
331  << Face_index_at_boundary[i][e] << std::endl;
332  }
333  }
334  }
335 
336 
337  // Lookup scheme has now been setup
339 
340 
341 }
342 
343 }
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
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
void setup_boundary_element_info()
Definition: quad_mesh.h:86
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
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806