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