hermite_element_quad_mesh.template.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 // the mesh assembly functions for the hermite element quad mesh
31 #ifndef OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC
32 #define OOMPH_HERMITE_QUAD_MESH_TEMPLATE_CC
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 //oomph-lib includes
42 
43 
44 namespace oomph
45 {
46 
47 
48 //=============================================================================
49 /// \short sets the generalised position of the node (i.e. - x_i, dx_i/ds_0,
50 /// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the
51 /// node from which its position can be determined.
52 //=============================================================================
53 template<class ELEMENT>
54 void HermiteQuadMesh<ELEMENT>::set_position_of_node(const unsigned& node_num_x,
55  const unsigned& node_num_y,
56  Node* node_pt)
57 {
58 
59  // get the generalised macro element position of the node
60  DenseMatrix<double> m_gen(4,2);
61  generalised_macro_element_position_of_node(node_num_x,node_num_y,m_gen);
62 
63  // copy the macro element coordinates
64  Vector<double> node_macro_position(2);
65  node_macro_position[0] = m_gen(0,0);
66  node_macro_position[1] = m_gen(0,1);
67 
68  // get the global position of the nodes
69  Vector<double> x_node(2);
70  Domain_pt->macro_element_pt(0)->macro_map(node_macro_position,x_node);
71 
72  // get the jacobian of the mapping from macro to eulerian coordinates
73  DenseMatrix<double> jacobian_MX(2,2);
74  Domain_pt->macro_element_pt(0)->
75  assemble_macro_to_eulerian_jacobian(node_macro_position,jacobian_MX);
76 
77 
78  // get the jacobian2 of the mapping from macro to eulerian coordinates
79  DenseMatrix<double> jacobian2_MX(3,2);
80  Domain_pt->macro_element_pt(0)->
81  assemble_macro_to_eulerian_jacobian2(node_macro_position,jacobian2_MX);
82 
83  // set x_0
84  node_pt->x_gen(0,0) = x_node[0];
85  // set x_1
86  node_pt->x_gen(0,1) = x_node[1];
87 
88  // set dxi/dsj for i = 0,1 and j = 0,1
89  // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
90  for (unsigned i = 0; i < 2; i++)
91  {
92  for (unsigned j = 0; j < 2; j++)
93  {
94  node_pt->x_gen(j+1,i) = m_gen(j+1,0)*jacobian_MX(0,i)
95  + m_gen(j+1,1)*jacobian_MX(1,i);
96  }
97  }
98 
99  // set d2xi/ds0ds1 for i = 0,1
100  // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
101  // + d2m1/ds0ds1*dxi/dm1
102  // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
103  // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
104  for (unsigned i = 0; i < 2; i++)
105  {
106  node_pt->x_gen(3,i) = m_gen(3,0)*jacobian_MX(0,i)
107  + m_gen(3,1)*jacobian_MX(1,i)
108  + m_gen(1,0)*(m_gen(2,0)*jacobian2_MX(0,i) +
109  m_gen(2,1)*jacobian2_MX(2,i))
110  + m_gen(1,1)*(m_gen(2,0)*jacobian2_MX(2,i) +
111  m_gen(2,1)*jacobian2_MX(1,i));
112  }
113 }
114 
115 
116  //=============================================================================
117 /// \short sets the generalised position of the node (i.e. - x_i, dx_i/ds_0,
118 /// dx_i/ds_1 & d2x_i/ds_0ds_1 for i = 1,2). Takes the x,y coordinates of the
119 /// node from which its position can be determined. Also sets coordinates
120 /// on boundary vector for the node to be the generalised position of the node
121 /// in macro element coordinates
122 //=============================================================================
123 template<class ELEMENT>
125 set_position_of_boundary_node(const unsigned& node_num_x,
126  const unsigned& node_num_y,
127  BoundaryNode<Node>* node_pt)
128 {
129 
130  // get the generalised macro element position of the node
131  DenseMatrix<double> m_gen(4,2,0.0);
132  generalised_macro_element_position_of_node(node_num_x,node_num_y,m_gen);
133 
134  // copy the macro element coordinates
135  Vector<double> node_macro_position(2);
136  node_macro_position[0] = m_gen(0,0);
137  node_macro_position[1] = m_gen(0,1);
138 
139  // get the global position of the nodes
140  Vector<double> x_node(2);
141  Domain_pt->macro_element_pt(0)->macro_map(node_macro_position,x_node);
142 
143  // get the jacobian of the mapping from macro to eulerian coordinates
144  DenseMatrix<double> jacobian_MX(2,2);
145  Domain_pt->macro_element_pt(0)->
146  assemble_macro_to_eulerian_jacobian(node_macro_position,jacobian_MX);
147 
148  // get the jacobian2 of the mapping from macro to eulerian coordinates
149  DenseMatrix<double> jacobian2_MX(3,2);
150  Domain_pt->macro_element_pt(0)->
151  assemble_macro_to_eulerian_jacobian2(node_macro_position,jacobian2_MX);
152 
153  // set x_0
154  node_pt->x_gen(0,0) = x_node[0];
155  // set x_1
156  node_pt->x_gen(0,1) = x_node[1];
157 
158  // set dxi/dsj for i = 0,1 and j = 0,1
159  // computation : dxi/dsj = dm0/dsj*dxi/dm0 + dm1/dsj*dxi/dm1
160  for (unsigned i = 0; i < 2; i++)
161  {
162  for (unsigned j = 0; j < 2; j++)
163  {
164  node_pt->x_gen(j+1,i) = m_gen(j+1,0)*jacobian_MX(0,i)
165  + m_gen(j+1,1)*jacobian_MX(1,i);
166  }
167  }
168 
169  // set d2xi/ds0ds1 for i = 0,1
170  // d2xi/ds0ds1 = d2m0/ds0ds1*dxi/dm0
171  // + d2m1/ds0ds1*dxi/dm1
172  // + (dm0/ds1*d2xi/dm0^2 + dm1/ds1*d2xi/dm0dm1)*dm0/ds0
173  // + (dm0/ds1*d2xi/dm0dm1 + dm1/ds1*d2xi/dm1^2)*dm1/ds0
174  for (unsigned i = 0; i < 2; i++)
175  {
176  node_pt->x_gen(3,i) = m_gen(3,0)*jacobian_MX(0,i)
177  + m_gen(3,1)*jacobian_MX(1,i)
178  + m_gen(1,0)*(m_gen(2,0)*jacobian2_MX(0,i) +
179  m_gen(2,1)*jacobian2_MX(2,i))
180  + m_gen(1,1)*(m_gen(2,0)*jacobian2_MX(2,i) +
181  m_gen(2,1)*jacobian2_MX(1,i));
182  }
183 
184 
185  // pass generalised macro element position to vector to pass to boundary
186  // node
187  for (unsigned b = 0; b < 4; b++)
188  {
189  if (node_pt->is_on_boundary(b))
190  {
191  Vector<double> boundary_macro_position(2,0);
192  boundary_macro_position[0] = m_gen(0,b%2);
193  boundary_macro_position[1] = m_gen(1+b%2,b%2);
194  node_pt->set_coordinates_on_boundary(b,boundary_macro_position);
195  }
196  }
197 }
198 
199 
200 //=============================================================================
201 /// \short computes the generalised position of the node at position
202 /// (node_num_x, node_num_y) in the macro element coordinate scheme.
203 /// index 0 of m_gen : 0 - m_i
204 /// 1 - dm_i/ds_0
205 /// 2 - dm_i/ds_1
206 /// 3 - d2m_i/ds_0ds_1 (where i is index 1 of m_gen)
207 //=============================================================================
208 template<class ELEMENT>
210 (const unsigned& node_num_x,const unsigned& node_num_y,
211  DenseMatrix<double>& m_gen)
212 {
213  // obtain position of node
214  Vector<double> s_macro_N(2);
215  macro_coordinate_position(node_num_x,node_num_y,s_macro_N);
216 
217  // set position of node in macro element coordinates
218  m_gen(0,0) = s_macro_N[0];
219  m_gen(0,1) = s_macro_N[1];
220 
221  // if on left hand edge
222  if (node_num_x == 0)
223  {
224 
225  // first dm0ds0 and dm1ds0
226 
227  // get position of right node
228  Vector<double> s_macro_R(2);
229  macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
230 
231  // compute dm0ds0
232  m_gen(1,0) = 0.5*(s_macro_R[0]-s_macro_N[0]);
233 
234  // compute dm1ds0
235  m_gen(1,1) = 0.5*(s_macro_R[1]-s_macro_N[1]);
236 
237  // now d2m0/ds0ds1 and d2m1/ds0ds1
238 
239  // if lower left hand corner
240  if (node_num_y == 0)
241  {
242 
243  // get position of upper right node
244  Vector<double> s_macro_UR(2);
245  macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
246 
247  // get position of upper node
248  Vector<double> s_macro_U(2);
249  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
250 
251  // compute d2m0/ds0ds1
252  m_gen(3,0) = 0.5 * (0.5*(s_macro_UR[0] - s_macro_U[0]) - m_gen(1,0));
253 
254  // compute d2m1/ds0ds1
255  m_gen(3,1) = 0.5 * (0.5*(s_macro_UR[1] - s_macro_U[1]) - m_gen(1,1));
256  }
257 
258  // else if upper left hand corner
259  else if (node_num_y == Nelement[1])
260  {
261 
262  // get position of lower right node
263  Vector<double> s_macro_DR(2);
264  macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
265 
266  // get position of lower node
267  Vector<double> s_macro_D(2);
268  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
269 
270  // compute d2m0/ds0ds1
271  m_gen(3,0) = 0.5 * (m_gen(1,0) - 0.5*(s_macro_DR[0]-s_macro_D[0]));
272 
273  // compute d2m0/ds0ds1
274  m_gen(3,1) = 0.5 * (m_gen(1,1) - 0.5*(s_macro_DR[1]-s_macro_D[1]));
275  }
276 
277  // else left hand edge (not corner)
278  else
279  {
280 
281  // get position of lower right node
282  Vector<double> s_macro_DR(2);
283  macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
284 
285  // get position of lower node
286  Vector<double> s_macro_D(2);
287  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
288 
289  // get position of upper right node
290  Vector<double> s_macro_UR(2);
291  macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
292 
293  // get position of upper node
294  Vector<double> s_macro_U(2);
295  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
296 
297  // compute d2m0/ds0ds1
298  m_gen(3,0) = 0.5*std::min(m_gen(1,0)-0.5*(s_macro_DR[0]-s_macro_D[0]),
299  0.5*(s_macro_UR[0]-s_macro_U[0])-m_gen(1,0));
300 
301  // compute d2m1/ds0ds1
302  m_gen(3,1) = 0.5*std::min(m_gen(1,1)-0.5*(s_macro_DR[1]-s_macro_D[1]),
303  0.5*(s_macro_UR[1]-s_macro_U[1])-m_gen(1,1));
304  }
305  }
306 
307  // if on right hand edge
308  else if (node_num_x == Nelement[0])
309  {
310 
311  // first dm0ds0 and dm1ds0
312 
313  // get position of left node
314  Vector<double> s_macro_L(2);
315  macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
316 
317  // compute dm0ds0
318  m_gen(1,0) = 0.5*(s_macro_N[0]-s_macro_L[0]);
319 
320  // compute dm1ds0
321  m_gen(1,1) = 0.5*(s_macro_N[1]-s_macro_L[1]);
322 
323  // now d2m0/ds0ds1 and d2m1/ds0ds1
324 
325  // if lower right hand corner
326  if (node_num_y == 0)
327  {
328 
329  // get position of upper left node
330  Vector<double> s_macro_UL(2);
331  macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);
332 
333  // get position of upper node
334  Vector<double> s_macro_U(2);
335  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
336 
337  // compute d2m0/ds0ds1
338  m_gen(3,0) = 0.5 * (0.5*(s_macro_U[0] - s_macro_UL[0]) - m_gen(1,0));
339 
340  // compute d2m1/ds0ds1
341  m_gen(3,1) = 0.5 * (0.5*(s_macro_U[1] - s_macro_UL[1]) - m_gen(1,1));
342  }
343 
344  // else if upper right hand corner
345  else if (node_num_y == Nelement[1])
346  {
347 
348  // get position of lower left node
349  Vector<double> s_macro_DL(2);
350  macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
351 
352  // get position of lower node
353  Vector<double> s_macro_D(2);
354  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
355 
356  // compute d2m0/ds0ds1
357  m_gen(3,0) = 0.5 * (m_gen(1,0) - 0.5*(s_macro_D[0]-s_macro_DL[0]));
358 
359  // compute d2m0/ds0ds1
360  m_gen(3,1) = 0.5 * (m_gen(1,1) - 0.5*(s_macro_D[1]-s_macro_DL[1]));
361  }
362 
363  // else left hand edge (not corner)
364  else
365  {
366 
367 
368  // get position of lower left node
369  Vector<double> s_macro_DL(2);
370  macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
371 
372  // get position of lower node
373  Vector<double> s_macro_D(2);
374  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
375 
376  // get position of upper left node
377  Vector<double> s_macro_UL(2);
378  macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);
379 
380  // get position of upper node
381  Vector<double> s_macro_U(2);
382  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
383 
384  // compute d2m0/ds0ds1
385  m_gen(3,0) = 0.5*std::min(m_gen(1,0)-0.5*(s_macro_D[0]-s_macro_DL[0]),
386  0.5*(s_macro_U[0]-s_macro_UL[0])-m_gen(1,0));
387 
388  // compute d2m0/ds0ds1
389  m_gen(3,1) = 0.5*std::min(m_gen(1,1)-0.5*(s_macro_D[1]-s_macro_DL[1]),
390  0.5*(s_macro_U[1]-s_macro_UL[1])-m_gen(1,1));
391  }
392  }
393  //
394  else
395  {
396  // get position of left node
397  Vector<double> s_macro_L(2);
398  macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
399 
400  // get position of right node
401  Vector<double> s_macro_R(2);
402  macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
403 
404  // compute dm0/ds0
405  m_gen(1,0) = 0.5*std::min(s_macro_N[0]-s_macro_L[0],
406  s_macro_R[0]-s_macro_N[0]);
407 
408  //compute dm1/ds0
409  Vector<double> node_space(2);
410  node_space[0] = s_macro_N[1]-s_macro_L[1];
411  node_space[1] = s_macro_R[1]-s_macro_N[1];
412  // set nodal dof
413  if (node_space[0] > 0 && node_space[1] > 0)
414  m_gen(1,1) = 0.5 * std::min(node_space[0],node_space[1]);
415  else if (node_space[0] < 0 && node_space[1] < 0)
416  m_gen(1,1) = 0.5 * std::max(node_space[0],node_space[1]);
417  else
418  m_gen(1,1) = 0;
419  }
420 
421  // if node in lower row
422  if (node_num_y == 0)
423  {
424 
425  // now dm0ds1 and dm1ds1
426 
427  // get position of upper node
428  Vector<double> s_macro_U(2);
429  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
430 
431  // compute dm0ds1
432  m_gen(2,0) = 0.5*(s_macro_U[0]-s_macro_N[0]);
433 
434  // compute dm1ds1
435  m_gen(2,1) = 0.5*(s_macro_U[1]-s_macro_N[1]);
436 
437  // and if not corner node
438  if (node_num_x > 0 && node_num_x < Nelement[0])
439  {
440 
441  // now dm0/ds0ds1 and dm1/ds0ds1
442 
443  // get position of upper left node
444  Vector<double> s_macro_UL(2);
445  macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);
446 
447  // get position of left node
448  Vector<double> s_macro_L(2);
449  macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
450 
451  // get position of right node
452  Vector<double> s_macro_R(2);
453  macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
454 
455  // get position of upper right node
456  Vector<double> s_macro_UR(2);
457  macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
458 
459  // set dm0/ds0ds1
460  m_gen(3,0) = 0.5*std::min(m_gen(2,0)-0.5*(s_macro_UL[0]-s_macro_L[0]),
461  0.5*(s_macro_UR[0]-s_macro_R[0])-m_gen(2,0));
462 
463  // set dm1/ds0ds1
464  m_gen(3,1) = 0.5*std::min(m_gen(2,1)-0.5*(s_macro_UL[1]-s_macro_L[1]),
465  0.5*(s_macro_UR[1]-s_macro_R[1])-m_gen(2,1));
466 
467  }
468  }
469 
470  // else if node in upper row
471  else if (node_num_y == Nelement[1])
472  {
473 
474  // now dm0ds1 and dm1ds1
475 
476  // get position of lower node
477  Vector<double> s_macro_D(2);
478  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
479 
480  // compute dm0ds1
481  m_gen(2,0) = 0.5*(s_macro_N[0]-s_macro_D[0]);
482 
483  // compute dm1ds1
484  m_gen(2,1) = 0.5*(s_macro_N[1]-s_macro_D[1]);
485 
486  // and if not corner node
487  if (node_num_x > 0 && node_num_x < Nelement[0])
488  {
489 
490  // now dm0/ds0ds1 and dm1/ds0ds1
491 
492  // get position of lower left node
493  Vector<double> s_macro_DL(2);
494  macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
495 
496  // get position of left node
497  Vector<double> s_macro_L(2);
498  macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
499 
500  // get position of right node
501  Vector<double> s_macro_R(2);
502  macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
503 
504  // get position of lower right node
505  Vector<double> s_macro_DR(2);
506  macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
507 
508  // set dm0/ds0ds1
509  m_gen(3,0) = 0.5*std::min(m_gen(2,0)-0.5*(s_macro_L[0]-s_macro_DL[0]),
510  0.5*(s_macro_R[0]-s_macro_DR[0])-m_gen(2,0));
511 
512  // set dm1/ds0ds1
513  m_gen(3,1) = 0.5*std::min(m_gen(2,1)-0.5*(s_macro_L[1]-s_macro_DL[1]),
514  0.5*(s_macro_R[1]-s_macro_DR[1])-m_gen(2,1));
515  }
516  }
517  else
518  {
519 
520  // get position of upper node
521  Vector<double> s_macro_U(2);
522  macro_coordinate_position(node_num_x,node_num_y+1,s_macro_U);
523 
524  // get position of lower node
525  Vector<double> s_macro_D(2);
526  macro_coordinate_position(node_num_x,node_num_y-1,s_macro_D);
527 
528  //compute dm0/ds1
529  Vector<double> node_space(2);
530  node_space[0] = s_macro_N[0]-s_macro_D[0];
531  node_space[1] = s_macro_U[0]-s_macro_N[0];
532  // set nodal coordinate
533  if (node_space[0] > 0 && node_space[1] > 0)
534  m_gen(2,0) = 0.5 * std::min(node_space[0],node_space[1]);
535  else if (node_space[0] < 0 && node_space[1] < 0)
536  m_gen(2,0) = 0.5 * std::max(node_space[0],node_space[1]);
537  else
538  m_gen(2,0) = 0;
539 
540  // compute dm1/ds1
541  m_gen(2,1) = 0.5*std::min(s_macro_N[1]-s_macro_D[1],
542  s_macro_U[1]-s_macro_N[1]);
543 
544  // for interior nodes
545  if (node_num_x > 0 && node_num_x < Nelement[0])
546  {
547 
548  // get position of left node
549  Vector<double> s_macro_L(2);
550  macro_coordinate_position(node_num_x-1,node_num_y,s_macro_L);
551 
552  // get position of upper left node
553  Vector<double> s_macro_UL(2);
554  macro_coordinate_position(node_num_x-1,node_num_y+1,s_macro_UL);
555 
556  // get position of upper right node
557  Vector<double> s_macro_UR(2);
558  macro_coordinate_position(node_num_x+1,node_num_y+1,s_macro_UR);
559 
560  // get position of right node
561  Vector<double> s_macro_R(2);
562  macro_coordinate_position(node_num_x+1,node_num_y,s_macro_R);
563 
564  // get position of lower right node
565  Vector<double> s_macro_DR(2);
566  macro_coordinate_position(node_num_x+1,node_num_y-1,s_macro_DR);
567 
568  // get position of lower left node
569  Vector<double> s_macro_DL(2);
570  macro_coordinate_position(node_num_x-1,node_num_y-1,s_macro_DL);
571 
572  Vector<double> node_space(2);
573  // comute dm0/ds0ds1 wrt to node above and below
574  node_space[0] = m_gen(1,0) - 0.5*std::min(s_macro_D[0]-s_macro_DL[0],
575  s_macro_DR[0]-s_macro_D[0]);
576  node_space[1] = 0.5*std::min(s_macro_U[0]-s_macro_UL[0],
577  s_macro_UR[0]-s_macro_U[0]) - m_gen(1,0);
578  // set nodal dof
579  if (node_space[0] > 0 && node_space[1] > 0)
580  m_gen(3,0) = 0.5 * std::min(node_space[0],node_space[1]);
581  else if (node_space[0] < 0 && node_space[1] < 0)
582  m_gen(3,0) = 0.5 * std::max(node_space[0],node_space[1]);
583  else
584  m_gen(3,0) = 0;
585 
586  // comute dm1/ds0ds1 wrt node left and right
587  node_space[0] = m_gen(2,1) - 0.5*std::min(s_macro_L[0]-s_macro_DL[0],
588  s_macro_UL[0]-s_macro_L[0]);
589  node_space[1] = 0.5*std::min(s_macro_R[0]-s_macro_DR[0],
590  s_macro_UR[0]-s_macro_R[0]) - m_gen(2,1);
591  // set nodal dof
592  if (node_space[0] > 0 && node_space[1] > 0)
593  m_gen(3,1) = 0.5 * std::min(node_space[0],node_space[1]);
594  else if (node_space[0] < 0 && node_space[1] < 0)
595  m_gen(3,1) = 0.5 * std::max(node_space[0],node_space[1]);
596  else
597  m_gen(3,1) = 0;
598  }
599  }
600  }
601 
602 
603 //=============================================================================
604 /// \short Generic mesh construction function to build the mesh
605 //=============================================================================
606 template<class ELEMENT>
608 {
609  // Mesh can only be built with 2D QHermiteElements.
610  MeshChecker::assert_geometric_element<QHermiteElementBase,ELEMENT>(2);
611 
612  //Set the number of boundaries
613  set_nboundary(4);
614 
615  //Allocate the store for the elements
616  Element_pt.resize(Nelement[0]*Nelement[1]);
617 
618  //Allocate the memory for the first element
619  Element_pt[0] = new ELEMENT;
620  finite_element_pt(0)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
621 
622  //Can now allocate the store for the nodes
623  Node_pt.resize((1 + Nelement[0])*(1 + Nelement[1]));
624 
625  //Set up geometrical data
626  unsigned long node_count=0;
627 
628  //Now assign the topology
629  //Boundaries are numbered 0 1 2 3 from the bottom proceeding anticlockwise
630  //Pinned value are denoted by an integer value 1
631  //Thus if a node is on two boundaries, ORing the values of the
632  //boundary conditions will give the most restrictive case (pinning)
633 
634 
635  // we first create the lowest row of elements
636  // ##########################################
637  // ##########################################
638 
639 
640  // FIRST ELEMENT - lower left hand corner
641  // **************************************
642 
643 
644  // LOWER LEFT HAND NODE
645 
646  //Set the corner node
647  //Allocate memory for the node
648  Node_pt[node_count] =
649  finite_element_pt(0)->construct_boundary_node(0,time_stepper_pt);
650 
651  //Push the node back onto boundaries
652  add_boundary_node(0,Node_pt[node_count]);
653  add_boundary_node(3,Node_pt[node_count]);
654 
655  // set the position of the boundary node
656  set_position_of_boundary_node
657  (0,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
658 
659  //Increment the node number
660  node_count++;
661 
662  // LOWER RIGHT HAND NODE
663 
664  //Allocate memory for the node
665  Node_pt[node_count] =
666  finite_element_pt(0)->construct_boundary_node(1,time_stepper_pt);
667 
668  //Push the node back onto boundaries
669  add_boundary_node(0,Node_pt[node_count]);
670 
671  //If we only have one column then the RHS node is on the right boundary
672  if(Nelement[0]==1)
673  {add_boundary_node(1,Node_pt[node_count]);}
674 
675  // set the position of the node
676  set_position_of_boundary_node
677  (1,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
678 
679  //Increment the node number
680  node_count++;
681 
682  // UPPER LEFT HAND NODE
683 
684  //Allocate memory for the node
685  Node_pt[node_count] =
686  finite_element_pt(0)->construct_boundary_node(2,time_stepper_pt);
687 
688  //Push the node back onto boundaries
689  add_boundary_node(3,Node_pt[node_count]);
690 
691  //If we only have one row, then the top node is on the top boundary
692  if(Nelement[1]==1)
693  {add_boundary_node(2,Node_pt[node_count]);}
694 
695  // set the position of the boundary node
696  set_position_of_boundary_node
697  (0,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
698 
699  //Increment the node number
700  node_count++;
701 
702  // UPPER RIGHT NODE
703 
704  //Allocate the memory for the node
705  Node_pt[node_count] =
706  finite_element_pt(0)->construct_node(3,time_stepper_pt);
707 
708  //If we only have one column then the RHS node is on the right boundary
709  if(Nelement[0]==1)
710  {add_boundary_node(1,Node_pt[node_count]);}
711  //If we only have one row, then the top node is on the top boundary
712  if(Nelement[1]==1)
713  {add_boundary_node(2,Node_pt[node_count]);}
714 
715  // if the node is a boundary node
716  if (Nelement[0]==1 || Nelement[1]==1)
717  {
718  // set the position of the boundary node
719  set_position_of_boundary_node
720  (1,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
721  }
722  else
723  {
724  // set the position of the node
725  set_position_of_node(1,1,Node_pt[node_count]);
726  }
727 
728  //Increment the node number
729  node_count++;
730 
731  // END OF FIRST ELEMENT
732 
733 
734  // CENTRE OF FIRST ROW OF ELEMENTS
735  // *******************************
736 
737 
738  //Now loop over the first row of elements, apart from final element
739  for(unsigned j=1;j<(Nelement[0]-1);j++)
740  {
741  //Allocate memory for new element
742  Element_pt[j] = new ELEMENT;
743  finite_element_pt(j)->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
744 
745  // LOWER LEFT NODE
746 
747  //lower left hand node column of nodes is same as lower right hand node
748  //in neighbouring element
749  finite_element_pt(j)->node_pt(0) =
750  finite_element_pt(j-1)->node_pt(1);
751 
752  // LOWER RIGHT NODE
753 
754  //Allocate memory for the nodes
755  Node_pt[node_count] =
756  finite_element_pt(j)->construct_boundary_node(1,time_stepper_pt);
757 
758  //Push the node back onto boundaries
759  add_boundary_node(0,Node_pt[node_count]);
760 
761  // set the position of the boundary node
762  set_position_of_boundary_node
763  (j+1,0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
764 
765  //Increment the node number
766  node_count++;
767 
768  // UPPER LEFT NODE
769 
770  //First column of nodes is same as neighbouring element
771  finite_element_pt(j)->node_pt(2) =
772  finite_element_pt(j-1)->node_pt(3);
773 
774  // UPPER RIGHT NODE
775 
776  //Allocate memory for the nodes
777  Node_pt[node_count] =
778  finite_element_pt(j)->construct_node(3,time_stepper_pt);
779 
780  //If we only have one row, then the top node is on the top boundary
781  if(Nelement[0]==1)
782  {add_boundary_node(2,Node_pt[node_count]);}
783 
784  // set the position of the (boundary) node
785  if (Nelement[0]==1)
786  set_position_of_boundary_node
787  (j+1,1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
788  else
789  set_position_of_node(j+1,1,Node_pt[node_count]);
790 
791  //Increment the node number
792  node_count++;
793  }
794 
795  // FINAL ELEMENT IN FIRST ROW (lower right corner element)
796  // **************************
797 
798 
799  //Only allocate if there is more than one element in the row
800  if(Nelement[0] > 1)
801  {
802  //Allocate memory for element
803  Element_pt[Nelement[0]-1] = new ELEMENT;
804  finite_element_pt(Nelement[0]-1)->
805  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
806 
807  // LOWER LEFT NODE
808 
809  //First column of nodes is same as neighbouring element
810  finite_element_pt(Nelement[0]-1)->node_pt(0) =
811  finite_element_pt(Nelement[0]-2)->node_pt(1);
812 
813  // LOWER RIGHT NODE
814 
815  //If we have an Xperiodic mesh then the final node is the first node in
816  //the row
817  if(Xperiodic==true)
818  {
819  //Note that this is periodic in the x-direction
820  finite_element_pt(Nelement[0]-1)->node_pt(1) =
821  finite_element_pt(0)->node_pt(0);
822  }
823  //Otherwise it's a new final node
824  else
825  {
826  //Allocate memory for the node
827  Node_pt[node_count] =
828  finite_element_pt(Nelement[0]-1)->
829  construct_boundary_node(1,time_stepper_pt);
830  }
831 
832  //Push the node back onto boundaries
833  add_boundary_node(0,Node_pt[node_count]);
834  add_boundary_node(1,Node_pt[node_count]);
835 
836  // set the position of the boundary node
837  set_position_of_boundary_node
838  (Nelement[0],0,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
839 
840  // if not periodic mesh
841  if(Xperiodic==false)
842  {
843  node_count++;
844  }
845 
846  // UPPER LEFT NODE
847 
848  // same as upper right node in element to left
849  finite_element_pt(Nelement[0]-1)->node_pt(2) =
850  finite_element_pt(Nelement[0]-2)->node_pt(3);
851 
852  //If we only have one row, then the top node is on the top boundary
853  if(Nelement[1]==1)
854  {add_boundary_node(2,Node_pt[node_count]);}
855 
856  // set the position of the (boundary) node
857  if (Nelement[1]==1)
858  set_position_of_boundary_node
859  (Nelement[0]-1,1,
860  dynamic_cast<BoundaryNode<Node> *>
861  (finite_element_pt(Nelement[0]-2)->node_pt(3)));
862 
863  // UPPER RIGHT NODE
864 
865  //If we have an Xperiodic mesh then the nodes in the final column are
866  //those in the first column
867  if(Xperiodic==true)
868  {
869  //Note that these are periodic in the x-direction
870  finite_element_pt(Nelement[0]-1)->node_pt(3) =
871  finite_element_pt(0)->node_pt(2);
872  }
873  //Otherwise we need new nodes for the final column
874  else
875  {
876  //Allocate memory for node
877  Node_pt[node_count] = finite_element_pt(Nelement[0]-1)->
878  construct_boundary_node(3,time_stepper_pt);
879  }
880 
881  //If we only have one row, then the top node is on the top boundary
882  if(Nelement[1]==1)
883  {add_boundary_node(2,Node_pt[node_count]);}
884 
885  // boundary 1 node
886  add_boundary_node(1,Node_pt[node_count]);
887 
888  // set the position of the boundary node
889  set_position_of_boundary_node
890  (Nelement[0],1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
891 
892  //Increment the node number
893  if(Xperiodic==false)
894  {
895  //Increment the node number
896  node_count++;
897  }
898  }
899 
900  // now create all remaining central rows
901  // #####################################
902  // #####################################
903 
904 
905  //Loop over remaining element rows in the mesh
906  for(unsigned i=1;i<(Nelement[1]-1);i++)
907  {
908 
909 
910  // set the first element in the row
911  // ********************************
912 
913 
914  //Allocate memory for element
915  Element_pt[Nelement[0]*i] = new ELEMENT;
916  finite_element_pt(Nelement[0]*i)
917  ->set_macro_elem_pt(Domain_pt->macro_element_pt(0));
918 
919  //The first row of nodes is copied from the element below
920  for(unsigned l=0;l<2;l++)
921  {
922  finite_element_pt(Nelement[0]*i)->node_pt(l) =
923  finite_element_pt(Nelement[0]*(i-1))->node_pt(2 + l);
924  }
925 
926  // UPPER LEFT HAND NODE
927 
928  //Allocate memory for node
929  Node_pt[node_count] = finite_element_pt(Nelement[0]*i)->
930  construct_boundary_node(2,time_stepper_pt);
931 
932  //Push the node back onto boundaries
933  add_boundary_node(3,Node_pt[node_count]);
934 
935  // set the position of the boundary node
936  set_position_of_boundary_node
937  (0,i+1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
938 
939  //Increment the node number
940  node_count++;
941 
942  // UPPER RIGHT HAND NODE
943 
944  //If we only have one column, the node could be on the boundary
945  if(Nelement[0]==1)
946  {
947  Node_pt[node_count] = finite_element_pt(Nelement[0]*i)->
948  construct_boundary_node(3,time_stepper_pt);
949  }
950  else
951  {
952  Node_pt[node_count] =
953  finite_element_pt(Nelement[0]*i)->
954  construct_node(3,time_stepper_pt);
955  }
956 
957  //If we only have one column then the RHS node is on the
958  //right boundary
959  if(Nelement[0]==1)
960  {add_boundary_node(1,Node_pt[node_count]);}
961 
962  // set the position of the (boundary) node
963  if (Nelement[0]==1)
964  set_position_of_boundary_node
965  (1,i+1,dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
966  else
967  set_position_of_node(1,i+1,Node_pt[node_count]);
968 
969  //Increment the node number
970  node_count++;
971 
972 
973  // loop over central elements in row
974  // *********************************
975 
976  //Now loop over the rest of the elements in the row, apart from the last
977  for(unsigned j=1;j<(Nelement[0]-1);j++)
978  {
979  //Allocate memory for new element
980  Element_pt[Nelement[0]*i+j] = new ELEMENT;
981  finite_element_pt(Nelement[0]*i+j)->
982  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
983 
984  // LOWER LEFT AND RIGHT NODES
985 
986  //The first row is copied from the lower element
987  for(unsigned l=0;l<2;l++)
988  {
989  finite_element_pt(Nelement[0]*i+j)->node_pt(l) =
990  finite_element_pt(Nelement[0]*(i-1)+j)->node_pt(2 + l);
991  }
992 
993  // UPPER LEFT NODE
994 
995  //First column is same as neighbouring element
996  finite_element_pt(Nelement[0]*i+j)->node_pt(2) =
997  finite_element_pt(Nelement[0]*i+(j-1))->node_pt(3);
998 
999  // UPPER RIGHT NODE
1000 
1001  //Allocate memory for the nodes
1002  Node_pt[node_count] =
1003  finite_element_pt(Nelement[0]*i+j)->
1004  construct_node(3,time_stepper_pt);
1005 
1006  // set position of node
1007  set_position_of_node(j+1,i+1,Node_pt[node_count]);
1008 
1009  //Increment the node number
1010  node_count++;
1011  }
1012 
1013 
1014  // final element in row
1015  // ********************
1016 
1017 
1018  //Only if there is more than one column
1019  if(Nelement[0] > 1)
1020  {
1021  //Allocate memory for element
1022  Element_pt[Nelement[0]*i+Nelement[0]-1] = new ELEMENT;
1023  finite_element_pt(Nelement[0]*i+Nelement[0]-1)->
1024  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1025 
1026  // LOWER LEFT AND RIGHT NODES
1027 
1028  //The first row is copied from the lower element
1029  for(unsigned l=0;l<2;l++)
1030  {
1031  finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(l) =
1032  finite_element_pt(Nelement[0]*(i-1)
1033  +Nelement[0]-1)->node_pt(2 + l);
1034  }
1035 
1036  // UPPER LEFT NODE
1037 
1038  // First node is same as neighbouring element
1039  finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(2) =
1040  finite_element_pt(Nelement[0]*i+Nelement[0]-2)->node_pt(3);
1041 
1042  // UPPER RIGHT NODE
1043 
1044  //If we have an Xperiodic mesh then this node is the same
1045  //as the first node
1046  if(Xperiodic==true)
1047  {
1048  //Allocate memory for node, periodic in x-direction
1049  finite_element_pt(Nelement[0]*i+Nelement[0]-1)->node_pt(3) =
1050  finite_element_pt(Nelement[0]*i)->node_pt(2);
1051  }
1052  //Otherwise allocate memory for a new node
1053  else
1054  {
1055  Node_pt[node_count] =
1056  finite_element_pt(Nelement[0]*i+Nelement[0]-1)->
1057  construct_boundary_node(3,time_stepper_pt);
1058  }
1059 
1060  //Push the node back onto boundaries
1061  add_boundary_node(1,Node_pt[node_count]);
1062 
1063  // set position of boundary node
1064  set_position_of_boundary_node
1065  (Nelement[0],i+1,
1066  dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
1067 
1068  //Increment the node number
1069  node_count++;
1070  }
1071  }
1072 
1073 
1074  // final element row
1075  // #################
1076  // #################
1077 
1078 
1079  // only if there is more than one row
1080  if(Nelement[1] > 1)
1081  {
1082 
1083  // first element in upper row (upper left corner)
1084  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1085 
1086  //Allocate memory for element
1087  Element_pt[Nelement[0]*(Nelement[1]-1)] = new ELEMENT;
1088  finite_element_pt(Nelement[0]*(Nelement[1]-1))->
1089  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1090 
1091  // LOWER LEFT AND LOWER RIGHT NODES
1092 
1093  //The first row of nodes is copied from the element below
1094  for(unsigned l=0;l<2;l++)
1095  {
1096  finite_element_pt(Nelement[0]*(Nelement[1]-1))->node_pt(l)
1097  = finite_element_pt(Nelement[0]*(Nelement[1]-2))->node_pt(2 + l);
1098  }
1099 
1100  // UPPER LEFT NODE
1101 
1102  //Allocate memory for node
1103  Node_pt[node_count] = finite_element_pt(Nelement[0]*(Nelement[1]-1))->
1104  construct_boundary_node(2,time_stepper_pt);
1105 
1106  //Push the node back onto boundaries
1107  add_boundary_node(2,Node_pt[node_count]);
1108  add_boundary_node(3,Node_pt[node_count]);
1109 
1110  // set position of boundary node
1111  set_position_of_boundary_node
1112  (0,Nelement[1],dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
1113 
1114  //Increment the node number
1115  node_count++;
1116 
1117  // UPPER RIGHT NODE
1118 
1119  //Allocate memory for the node
1120  Node_pt[node_count] =
1121  finite_element_pt(Nelement[0]*(Nelement[1]-1))->
1122  construct_boundary_node(3,time_stepper_pt);
1123 
1124  //Push the node back onto boundaries
1125  add_boundary_node(2,Node_pt[node_count]);
1126 
1127  // set position of boundary node
1128  set_position_of_boundary_node
1129  (1,Nelement[1],dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
1130 
1131  //Increment the node number
1132  node_count++;
1133 
1134 
1135  // loop over central elements of upper element row
1136  // ***********************************************
1137 
1138 
1139  //Now loop over the rest of the elements in the row, apart from the last
1140  for(unsigned j=1;j<(Nelement[0]-1);j++)
1141  {
1142 
1143  //Allocate memory for element
1144  Element_pt[Nelement[0]*(Nelement[1]-1)+j] = new ELEMENT;
1145  finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->
1146  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1147 
1148  // LOWER LEFT AND LOWER RIGHT NODES
1149 
1150  //The first row is copied from the lower element
1151  for(unsigned l=0;l<2;l++)
1152  {
1153  finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->node_pt(l) =
1154  finite_element_pt(Nelement[0]*(Nelement[1]-2)+j)->node_pt(2 + l);
1155  }
1156 
1157  // UPPER LEFT NODE
1158 
1159  //First column is same as neighbouring element
1160  finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->node_pt(2)
1161  = finite_element_pt(Nelement[0]*(Nelement[1]-1)+(j-1))->node_pt(3);
1162 
1163  // UPPER RIGHT NODE
1164 
1165  //Allocate memory for node
1166  Node_pt[node_count] =
1167  finite_element_pt(Nelement[0]*(Nelement[1]-1)+j)->
1168  construct_boundary_node(3,time_stepper_pt);
1169 
1170  //Push the node back onto boundaries
1171  add_boundary_node(2,Node_pt[node_count]);
1172 
1173  // set position of boundary node
1174  set_position_of_boundary_node
1175  (j+1,Nelement[1],
1176  dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
1177 
1178  //Increment the node number
1179  node_count++;
1180 
1181  } //End of loop over central elements in row
1182 
1183 
1184  // final element in row (upper right corner)
1185  // *****************************************
1186 
1187  //Only if there is more than one column
1188  if(Nelement[0] > 1)
1189  {
1190 
1191  //Allocate memory for element
1192  Element_pt[Nelement[0]*(Nelement[1]-1)+Nelement[0]-1] = new ELEMENT;
1193  finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)->
1194  set_macro_elem_pt(Domain_pt->macro_element_pt(0));
1195 
1196  // LOWER LEFT AND LOWER RIGHT NODES
1197 
1198  //The first row is copied from the lower element
1199  for(unsigned l=0;l<2;l++)
1200  {
1201  finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)->
1202  node_pt(l) = finite_element_pt(Nelement[0]*(Nelement[1]-2)+
1203  Nelement[0]-1)->node_pt(2 + l);
1204  }
1205 
1206  // UPPER LEFT NODE
1207 
1208  //First column is same as neighbouring element
1209  finite_element_pt(Nelement[0]*(Nelement[1]-1)+Nelement[0]-1)
1210  ->node_pt(2) = finite_element_pt(Nelement[0]*(Nelement[1]-1)+
1211  Nelement[0]-2)->node_pt(3);
1212 
1213  // UPPER RIGHT NODE
1214 
1215  //If we have an Xperiodic mesh, the node must be copied from
1216  //the first column
1217  if(Xperiodic==true)
1218  {
1219  //Allocate memory for node, periodic in x-direction
1220  finite_element_pt(Nelement[0]*(Nelement[1]-1)
1221  +Nelement[0]-1)->node_pt(3) =
1222  finite_element_pt(Nelement[0]*(Nelement[1]-1))->node_pt(2);
1223  }
1224  //Otherwise, allocate new memory for node
1225  else
1226  {
1227  Node_pt[node_count] = finite_element_pt(Nelement[0]*(Nelement[1]-1)+
1228  Nelement[0]-1)->
1229  construct_boundary_node(3,time_stepper_pt);
1230  }
1231 
1232  //Push the node back onto boundaries
1233  add_boundary_node(1,Node_pt[node_count]);
1234  add_boundary_node(2,Node_pt[node_count]);
1235 
1236  // set position of boundary node
1237  set_position_of_boundary_node
1238  (Nelement[0],Nelement[1],
1239  dynamic_cast<BoundaryNode<Node> *>(Node_pt[node_count]));
1240 
1241  //Increment the node number
1242  node_count++;
1243  }
1244  }
1245 
1246  // Setup boundary element lookup schemes
1247  setup_boundary_element_info();
1248 }
1249 
1250 
1251 
1252 //=============================================================================
1253 /// \short Setup lookup schemes which establish which elements are located next
1254 /// to which boundaries (Doc to outfile if it's open).
1255 /// Specific version for HermiteQuadMesh to ensure that the order of the
1256 /// elements in Boundary_element_pt matches the actual order along the
1257 /// boundary. This is required when hijacking the BiharmonicElement to apply
1258 /// the BiharmonicFluidBoundaryElement in
1259 /// BiharmonicFluidProblem::impose_traction_free_edge(...)
1260 //================================================================
1261 template <class ELEMENT>
1263 setup_boundary_element_info(std::ostream &outfile)
1264 {
1265 
1266  bool doc=false;
1267  if (outfile) doc=true;
1268 
1269  // Number of boundaries
1270  unsigned nbound=nboundary();
1271 
1272  // Wipe/allocate storage for arrays
1273  Boundary_element_pt.clear();
1274  Face_index_at_boundary.clear();
1275  Boundary_element_pt.resize(nbound);
1276  Face_index_at_boundary.resize(nbound);
1277 
1278  // Temporary vector of sets of pointers to elements on the boundaries:
1279  Vector<Vector<FiniteElement*> > vector_of_boundary_element_pt;
1280  vector_of_boundary_element_pt.resize(nbound);
1281 
1282  // Matrix map for working out the fixed local coord for elements on boundary
1284  boundary_identifier;
1285 
1286 
1287  // Loop over elements
1288  //-------------------
1289  unsigned nel=nelement();
1290  for (unsigned e=0;e<nel;e++)
1291  {
1292  // Get pointer to element
1293  FiniteElement* fe_pt=finite_element_pt(e);
1294 
1295  if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
1296 
1297  // Only include 2D elements! Some meshes contain interface elements too.
1298  if (fe_pt->dim()==2)
1299  {
1300  // Loop over the element's nodes and find out which boundaries they're on
1301  // ----------------------------------------------------------------------
1302  unsigned nnode_1d=fe_pt->nnode_1d();
1303 
1304  // Loop over nodes in order
1305  for (unsigned i0=0;i0<nnode_1d;i0++)
1306  {
1307  for (unsigned i1=0;i1<nnode_1d;i1++)
1308  {
1309  // Local node number
1310  unsigned j=i0+i1*nnode_1d;
1311 
1312  // Get pointer to vector of boundaries that this
1313  // node lives on
1314  std::set<unsigned>* boundaries_pt=0;
1315  fe_pt->node_pt(j)->get_boundaries_pt(boundaries_pt);
1316 
1317  // If the node lives on some boundaries....
1318  if (boundaries_pt!=0)
1319  {
1320  // Loop over boundaries
1321  //unsigned nbound=(*boundaries_pt).size();
1322  for (std::set<unsigned>::iterator it=boundaries_pt->begin();
1323  it != boundaries_pt->end();++it)
1324  {
1325  // Add pointer to finite element to set for the appropriate
1326  // boundary -- storage in set makes sure we don't count elements
1327  // multiple times
1328  unsigned temp_size = vector_of_boundary_element_pt[*it].size();
1329  bool temp_flag = true;
1330  for (unsigned temp_i = 0; temp_i < temp_size; temp_i++)
1331  {
1332  if (vector_of_boundary_element_pt[*it][temp_i] == fe_pt)
1333  temp_flag = false;
1334  }
1335  if (temp_flag)
1336  vector_of_boundary_element_pt[*it].push_back(fe_pt);
1337 
1338  // For the current element/boundary combination, create
1339  // a vector that stores an indicator which element boundaries
1340  // the node is located (boundary_identifier=-/+1 for nodes
1341  // on the left/right boundary; boundary_identifier=-/+2 for nodes
1342  // on the lower/upper boundary. We determine these indices
1343  // for all corner nodes of the element and add them to a vector
1344  // to a vector. This allows us to decide which face of the element
1345  // coincides with the boundary since the (quad!) element must
1346  // have exactly two corner nodes on the boundary.
1347  if (boundary_identifier(*it,fe_pt)==0)
1348  {
1349  boundary_identifier(*it,fe_pt)=new Vector<int>;
1350  }
1351 
1352  // Are we at a corner node?
1353  if (((i0==0)||(i0==nnode_1d-1))&&((i1==0)||(i1==nnode_1d-1)))
1354  {
1355  // Create index to represent position relative to s_0
1356  (*boundary_identifier(*it,fe_pt)).
1357  push_back(1*(2*i0/(nnode_1d-1)-1));
1358 
1359  // Create index to represent position relative to s_1
1360  (*boundary_identifier(*it,fe_pt)).
1361  push_back(2*(2*i1/(nnode_1d-1)-1));
1362  }
1363 
1364  }
1365  }
1366  //else
1367  // {
1368  // oomph_info << "...does not live on any boundaries " << std::endl;
1369  // }
1370 
1371  }
1372  }
1373  }
1374  //else
1375  //{
1376  //oomph_info << "Element " << e << " does not qualify" << std::endl;
1377  //}
1378  }
1379 
1380 
1381  // Now copy everything across into permanent arrays
1382  //-------------------------------------------------
1383 
1384  // Note: vector_of_boundary_element_pt contains all elements
1385  // that have (at least) one corner node on a boundary -- can't copy
1386  // them across into Boundary_element_pt
1387  // yet because some of them might have only one node on the
1388  // the boundary in which case they don't qualify as
1389  // boundary elements!
1390 
1391  // Loop over boundaries
1392  //---------------------
1393  for (unsigned i=0;i<nbound;i++)
1394  {
1395  // Number of elements on this boundary (currently stored in a set)
1396  unsigned nel=vector_of_boundary_element_pt[i].size();
1397 
1398  // Allocate storage for the coordinate identifiers
1399  Face_index_at_boundary[i].resize(nel);
1400 
1401  // Loop over elements that have at least one corner node on this boundary
1402  //-----------------------------------------------------------------------
1403  unsigned e_count=0;
1405  for (IT it=vector_of_boundary_element_pt[i].begin();
1406  it!=vector_of_boundary_element_pt[i].end();
1407  it++)
1408  {
1409  // Recover pointer to element
1410  FiniteElement* fe_pt=*it;
1411 
1412  // Initialise count for boundary identiers (-2,-1,1,2)
1413  std::map<int,unsigned> count;
1414 
1415  // Loop over coordinates
1416  for (unsigned ii=0;ii<2;ii++)
1417  {
1418  // Loop over upper/lower end of coordinates
1419  for (int sign=-1;sign<3;sign+=2)
1420  {
1421  count[(ii+1)*sign]=0;
1422  }
1423  }
1424 
1425  // Loop over boundary indicators for this element/boundary
1426  unsigned n_indicators=(*boundary_identifier(i,fe_pt)).size();
1427  for (unsigned j=0;j<n_indicators;j++)
1428  {
1429  count[(*boundary_identifier(i,fe_pt))[j] ]++;
1430  }
1431  delete boundary_identifier(i,fe_pt);
1432 
1433  // Determine the correct boundary indicator by checking that it
1434  // occurs twice (since two corner nodes of the element's boundary
1435  // need to be located on the domain boundary
1436  int indicator=-10;
1437 
1438  //Check that we're finding exactly one boundary indicator
1439  bool found=false;
1440 
1441  // Loop over coordinates
1442  for (unsigned ii=0;ii<2;ii++)
1443  {
1444  // Loop over upper/lower end of coordinates
1445  for (int sign=-1;sign<3;sign+=2)
1446  {
1447  if (count[(ii+1)*sign]==2)
1448  {
1449  // Check that we haven't found multiple boundaries
1450  if (found)
1451  {
1452  std::string error_message=
1453  "Trouble: Multiple boundary identifiers!\n";
1454  error_message +=
1455  "Elements should only have at most 2 corner ";
1456  error_message +=
1457  "nodes on any one boundary.\n";
1458 
1459  throw OomphLibError(
1460  error_message,
1461  OOMPH_CURRENT_FUNCTION,
1462  OOMPH_EXCEPTION_LOCATION);
1463  }
1464  found=true;
1465  indicator=(ii+1)*sign;
1466  }
1467  }
1468  }
1469 
1470  // Element has exactly two corner nodes on boundary
1471  if (found)
1472  {
1473  // Add to permanent storage
1474  Boundary_element_pt[i].push_back(fe_pt);
1475 
1476  // Now convert boundary indicator into information required
1477  // for FaceElements
1478  switch (indicator)
1479  {
1480 
1481  case -2:
1482 
1483  // s_1 is fixed at -1.0:
1484  Face_index_at_boundary[i][e_count]= -2;
1485  break;
1486 
1487  case -1:
1488 
1489  // s_0 is fixed at -1.0:
1490  Face_index_at_boundary[i][e_count] = -1;
1491  break;
1492 
1493 
1494  case 1:
1495 
1496  // s_0 is fixed at 1.0:
1497  Face_index_at_boundary[i][e_count] = 1;
1498  break;
1499 
1500  case 2:
1501 
1502  // s_1 is fixed at 1.0:
1503  Face_index_at_boundary[i][e_count] = 2;
1504  break;
1505 
1506  default:
1507 
1508  throw OomphLibError("Never get here",
1509  OOMPH_CURRENT_FUNCTION,
1510  OOMPH_EXCEPTION_LOCATION);
1511  }
1512 
1513  // Increment counter
1514  e_count++;
1515  }
1516 
1517  }
1518  }
1519 
1520 
1521 
1522  // Doc?
1523  //-----
1524  if (doc)
1525  {
1526  // Loop over boundaries
1527  for (unsigned i=0;i<nbound;i++)
1528  {
1529  unsigned nel=Boundary_element_pt[i].size();
1530  outfile << "Boundary: " << i
1531  << " is adjacent to " << nel
1532  << " elements" << std::endl;
1533 
1534  // Loop over elements on given boundary
1535  for (unsigned e=0;e<nel;e++)
1536  {
1537  FiniteElement* fe_pt=Boundary_element_pt[i][e];
1538  outfile << "Boundary element:" << fe_pt
1539  << " Face index of element along boundary is "
1540  << Face_index_at_boundary[i][e] << std::endl;
1541  }
1542  }
1543  }
1544 
1545 
1546  // Lookup scheme has now been setup yet
1547  Lookup_for_elements_next_boundary_is_setup=true;
1548 
1549 }
1550 }
1551 
1552 #endif
void set_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, Node *node_pt)
sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = ...
cstr elem_len * i
Definition: cfortran.h:607
bool is_on_boundary() const
Test whether the node lies on a boundary Final overload.
Definition: nodes.h:2361
virtual void build_mesh(TimeStepper *time_stepper_pt)
Generic mesh construction function to build the mesh.
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
e
Definition: cfortran.h:575
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
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
void set_position_of_boundary_node(const unsigned &node_num_x, const unsigned &node_num_y, BoundaryNode< Node > *node_pt)
sets the generalised position of the node (i.e. - x_i, dx_i/ds_0, dx_i/ds_1 & d2x_i/ds_0ds_1 for i = ...
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh...
Definition: nodes.h:67
void set_coordinates_on_boundary(const unsigned &b, const Vector< double > &boundary_zeta)
Set the vector of coordinates on mesh boundary b Final overload.
Definition: nodes.h:2392
double & x_gen(const unsigned &k, const unsigned &i)
Reference to the generalised position x(k,i).
Definition: nodes.h:1055
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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 generalised_macro_element_position_of_node(const unsigned &node_num_x, const unsigned &node_num_y, DenseMatrix< double > &m_gen)
computes the generalised position of the node at position (node_num_x, node_num_y) in the macro eleme...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219