extruded_cube_mesh_from_quad_mesh_with_macro_elements.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-2022 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// Oomph-lib headers
28
29/// //////////////////////////////////////////////////////////////////////
30/// //////////////////////////////////////////////////////////////////////
31/// //////////////////////////////////////////////////////////////////////
32
33namespace oomph
34{
35 //========================================================================
36 /// Get all the boundary information of an element using the
37 /// input (quad_mesh_pt) mesh. If the element lies on a boundary then
38 /// the user will be given the corresponding boundary index and the
39 /// index of the face of quad_el_pt attached to the boundary. If the
40 /// element does NOT lie on any boundaries, this function simply
41 /// returns a vector of size zero.
42 //========================================================================
43 template<class ELEMENT>
44 Vector<std::pair<unsigned, int>> ExtrudedCubeMeshFromQuadMesh<
45 ELEMENT>::get_element_boundary_information(QuadMeshBase* quad_mesh_pt,
46 FiniteElement* quad_el_pt)
47 {
48 // Allocate space for the boundary info (initialised to empty)
50
51 // Find the number of boundaries in the input mesh
52 unsigned n_boundary = quad_mesh_pt->nboundary();
53
54 // Loop over the boundaries of the mesh
55 for (unsigned b = 0; b < n_boundary; b++)
56 {
57 // How many elements are there on the b-th boundary?
58 unsigned n_boundary_element = quad_mesh_pt->nboundary_element(b);
59
60 // Loop over the elements on the b-th boundary
61 for (unsigned e = 0; e < n_boundary_element; e++)
62 {
63 // Is the e-th element on the b-th boundary the input element?
64 if (quad_el_pt == quad_mesh_pt->boundary_element_pt(b, e))
65 {
66 // Create a pair to hold the boundary index and element face index
67 std::pair<unsigned, int> boundary_and_face_index;
68
69 // Set the first entry (boundary index)
70 boundary_and_face_index.first = b;
71
72 // Set the second entry (face index)
73 boundary_and_face_index.second =
74 quad_mesh_pt->face_index_at_boundary(b, e);
75
76 // Add the boundary index to the boundary_info vector
77 boundary_info.push_back(boundary_and_face_index);
78 }
79 } // for (unsigned e=0;e<n_boundary_element;e++)
80 } // for (unsigned b=0;b<n_boundary;b++)
81
82 // Return the result
83 return boundary_info;
84 } // End of get_element_boundary_information
85
86 //========================================================================
87 /// Generic mesh construction. This function contains all the details of
88 /// the mesh generation process, including all the tedious loops, counting
89 /// spacing and boundary functions.
90 /// NOTE: The boundary number of the extruded mesh will follow the same
91 /// numbering as the input quad mesh. The newly created "front" and "back"
92 /// face of the 3D mesh will added after those boundaries. For example,
93 /// if the input mesh has 4 boundaries; b_0, b_1, b_2 & b_3 then the
94 /// extruded mesh will have 6 boundaries; eb_0, eb_1, eb_2, eb_3, eb_4 &
95 /// eb_5 where eb_4 is the "front" face and eb5 is the "back" face. The
96 /// boundaries eb_i here satisfy eb_i = (b_i) x [Zmin,Zmax].
97 //========================================================================
98 template<class ELEMENT>
100 QuadMeshBase* quad_mesh_pt, TimeStepper* time_stepper_pt)
101 {
102 // Mesh can only be built with 3D Qelements.
103 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
104
105 // Need at least 2 elements in the z-direction
106 if (Nz == 0)
107 {
108 // Create an ostringstream object to create an error message
109 std::ostringstream error_message;
110
111 // Create the error message
112 error_message << "Extrude2DQuadMeshTo3DCubeMesh needs at least two "
113 << "elements in each,\ncoordinate direction. You have "
114 << "specified Nz=" << Nz << std::endl;
115
116 // Throw an error
117 throw OomphLibError(
118 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
119 }
120
121 // Get the current (i.e. start) time
122 double t_mesh_setup_start = TimingHelpers::timer();
123
124 // Get the current (i.e. start) time
125 double t_mesh_construct_start = TimingHelpers::timer();
126
127 // Get the number of boundaries in the input mesh
128 unsigned n_boundary_in_quad_mesh = quad_mesh_pt->nboundary();
129
130 // The number of boundaries in the extruded mesh
131 unsigned n_boundary_in_extruded_mesh = n_boundary_in_quad_mesh + 2;
132
133 // Set the number of boundaries in the extruded mesh (just need
134 // to add the front and back face to the number of boundaries)
135 set_nboundary(n_boundary_in_extruded_mesh);
136
137 // Get the number of elements in the input mesh
138 unsigned n_element_in_quad_mesh = quad_mesh_pt->nelement();
139
140 // Allocate storage for the boundary information of each quad element.
141 // Takes up more memory but avoids recomputing it many, many times
142 Vector<Vector<std::pair<unsigned, int>>> boundary_information(
143 n_element_in_quad_mesh);
144
145 // Loop over all the elements in the quad mesh
146 for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
147 {
148 // Get a pointer to the j-th element in the mesh
149 FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(e);
150
151 // Get the boundary information
152 boundary_information[e] =
153 get_element_boundary_information(quad_mesh_pt, quad_el_pt);
154 }
155
156 // Allocate storage for the number of elements in the extruded mesh
157 Element_pt.resize(n_element_in_quad_mesh * Nz);
158
159 // Create the first element
160 ELEMENT* temp_el_pt = new ELEMENT;
161
162 // Read out the number of linear points in the element (in one direction)
163 unsigned n_node_1d = temp_el_pt->nnode_1d();
164
165 // Store the variable n_node_1d as the private variable, Np
166 N_node_1d = n_node_1d;
167
168 // Delete the element
169 delete temp_el_pt;
170
171 // Make it a null pointer
172 temp_el_pt = 0;
173
174 // Need the same number of nodes in one direction in the 2D and 3D element
175 if ((n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode_1d()) ||
176 (n_node_1d * n_node_1d != quad_mesh_pt->finite_element_pt(0)->nnode()))
177 {
178 // Create an ostringstream object to create an error message
179 std::ostringstream error_message;
180
181 // Create the error message
182 error_message
183 << "Extrude2DQuadMeshTo3DCubeMesh needs the number of "
184 << "nodes (in the 2D element) in one direction\n to match "
185 << "the number of nodes in one direction in the 3D element!";
186
187 // Throw an error
188 throw OomphLibError(
189 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
190 }
191
192 // Get the total number of nodes in the input mesh
193 unsigned long n_node_in_quad_mesh = quad_mesh_pt->nnode();
194
195 // Loop over the nodes in the spatial mesh
196 for (unsigned i = 0; i < n_node_in_quad_mesh; i++)
197 {
198 // Check to see if the i-th node is hanging
199 if (quad_mesh_pt->node_pt(i)->is_hanging())
200 {
201 // Create an output stream
202 std::ostringstream warning_message_stream;
203
204 // Create an error message
205 warning_message_stream << "Extrusion machinery is not currently able "
206 << "to deal with hanging nodes!" << std::endl;
207
208 // Throw an error
209 throw OomphLibError(warning_message_stream.str(),
210 OOMPH_CURRENT_FUNCTION,
211 OOMPH_EXCEPTION_LOCATION);
212 }
213 } // for (unsigned i=0;i<n_node_in_quad_mesh;i++)
214
215 // Number of nodes in the z-direction
216 unsigned n_node_in_z_direction = (1 + (n_node_1d - 1) * Nz);
217
218 // Get the total number of nodes in the input mesh
219 unsigned long n_node_in_cube_mesh =
220 n_node_in_quad_mesh * n_node_in_z_direction;
221
222 // Allocate storage for the number of nodes in the extruded mesh
223 Node_pt.resize(n_node_in_cube_mesh);
224
225 // Counter for the number of nodes in the extruded mesh
226 unsigned long node_count = 0;
227
228 // Loop over each element slice in the z-direction
229 for (unsigned i = 0; i < Nz; i++)
230 {
231 // Need a map at each element slice to tell us whether or not we've
232 // already created an edge (and thus the nodes on that edge). If the
233 // nodes have already been created then all we have to do is copy them
234 // over in reverse order
235 std::map<std::pair<Edge, double>, Vector<Node*>> edge_to_nodes_map;
236
237 // As several edges can attach to a single (corner) node we can
238 // accidentally create duplicate nodes (i.e. nodes that lie at
239 // the same spatial position) because the edges aren't covered in
240 // the required order (something we shouldn't rely on anyway!). To
241 // get around this we store the corner nodes of each quad element
242 // with its corresponding z-value. This will in turn return the
243 // corresponding node in the extruded mesh
244 std::map<std::pair<Node*, double>, Node*>
245 quad_corner_node_to_cube_node_map;
246
247 // Loop over all the elements in one slice
248 for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
249 {
250 // Get a pointer to the e-th element in the input mesh
251 FiniteElement* quad_el_pt = quad_mesh_pt->finite_element_pt(e);
252
253 // Get the index of the cube element in the extruded mesh
254 unsigned element_index = i * n_element_in_quad_mesh + e;
255
256 // Create the e-th element in the i-th slice of the extruded mesh
257 ELEMENT* cube_el_pt = new ELEMENT;
258
259 // Put it in global storage!
260 Element_pt[element_index] = cube_el_pt;
261
262 // Index of the local node number (in the quad element)
263 unsigned quad_local_node_number = 0;
264
265 // Index of the local node number (in the cube element)
266 unsigned cube_local_node_number = 0;
267
268 // Cache the boundary information associated with this element
269 Vector<std::pair<unsigned, int>> el_boundary_information =
270 boundary_information[e];
271
272 // Loop over the nodes in the third local coordinate direction
273 for (unsigned n_2 = 0; n_2 < n_node_1d; n_2++)
274 {
275 // Calculate the z-value at the n_2-th (z-)node in the i-th element
276 // slice
277 double z_value = z_spacing_function(i, n_2);
278
279 //---------------------------------------------------------------
280 // At a given node slice, there are three options:
281 // (1) If i>0 and n_2=0 then we already created the nodes
282 // (in the current slice) when we created the nodes in
283 // the previous slice;
284 // (2) If (1) doesn't hold we may have still have created
285 // certain nodes along a face/edge of an element that
286 // we previously created (in this node slice);
287 // (3) No nodes in this node slice have been created yet so
288 // set them all up!
289 // NOTE: if (1) occurs, we can skip to the next node slice but
290 // we can't if (2) occurs (or (3), obviously!).
291 //---------------------------------------------------------------
292 //---------------------------------------------------------------
293 // Case (1): If we're past the first element slice and we're on
294 // the first node slice (in the element) then copy the nodes from
295 // the appropriate element.
296 //---------------------------------------------------------------
297 if ((i > 0) && (n_2 == 0))
298 {
299 // Get the index of the cube element (to copy) in the extruded mesh
300 unsigned copy_element_index = (i - 1) * n_element_in_quad_mesh + e;
301
302 // Get a pointer to the element in the previous slice
303 // containing the node we want
304 FiniteElement* copy_cube_el_pt =
305 finite_element_pt(copy_element_index);
306
307 // Storage for the local node number in the element that we're
308 // copying
309 unsigned copy_cube_local_node_number = 0;
310
311 // Loop over the nodes in the second local coordinate direction
312 for (unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
313 {
314 // Loop over the nodes in the first local coordinate direction
315 for (unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
316 {
317 // Calculate the local node number (in the extruded mesh)
318 cube_local_node_number = n_0 + (n_node_1d * n_1);
319
320 // Calculate the local node number (in the input mesh)
321 quad_local_node_number = n_0 + (n_node_1d * n_1);
322
323 // Calculate the local node number (in the element we're
324 // copying it from)
325 copy_cube_local_node_number =
326 cube_local_node_number +
327 n_node_1d * n_node_1d * (n_node_1d - 1);
328
329 // Copy the node over
330 cube_el_pt->node_pt(cube_local_node_number) =
331 copy_cube_el_pt->node_pt(copy_cube_local_node_number);
332 } // for (unsigned n_0=0;n_0<n_node_1d;n_0++)
333 } // for (unsigned n_1=0;n_1<n_node_1d;n_1++)
334 }
335 //---------------------------------------------------------------
336 // Case (2) & (3): Loop over the edges in the current node slice
337 // and check if any of them have already been set up (and thus
338 // stored in edge_to_nodes_map). If they have then copy them over
339 // and construct all the other nodes. If they haven't then simply
340 // construct all of the nodes.
341 // NOTE: if an edge has been created already, we can't make the
342 // nodes first, assign the coordinates THEN check as it would
343 // be a wasteful operation. Instead, we get the nodes describing
344 // an edge from the input mesh (quad_mesh_pt) and pair it with
345 // the current z-value.
346 //---------------------------------------------------------------
347 else
348 {
349 // Create a vector of bools to tell us if we've set the nodes up
350 // in the current node slice
351 std::vector<bool> has_node_been_setup(n_node_1d * n_node_1d, false);
352
353 // Number of edges in a slice (N/E/S/W)
354 unsigned n_edge = 4;
355
356 // Storage for the edges
357 Vector<std::pair<Edge, double>> edges_and_z_value;
358
359 // Vector to tell us if we are copying any edges (to make sure
360 // we don't store the edge again later!)
361 std::vector<bool> has_edge_been_done(n_edge, false);
362
363 // What is the last edge? Convert this to an int so the compiler
364 // doesn't complain about comparisons between unsigned and signed
365 // integers...
366 int last_edge = int(QuadTreeNames::N + n_edge);
367
368 // Loop over the edges
369 for (int i_edge = QuadTreeNames::N; i_edge < last_edge; i_edge++)
370 {
371 // Pointer for the first node associated with this edge
372 Node* node1_pt = 0;
373
374 // Pointer for the second node associated with this edge
375 Node* node2_pt = 0;
376
377 // Find the corner nodes associated with this edge
378 switch (i_edge)
379 {
380 case QuadTreeNames::N:
381 // The first node
382 node1_pt = quad_el_pt->node_pt((n_node_1d * n_node_1d) - 1);
383
384 // The second node
385 node2_pt = quad_el_pt->node_pt(n_node_1d * (n_node_1d - 1));
386
387 // Break
388 break;
389 case QuadTreeNames::E:
390 // The first node
391 node1_pt = quad_el_pt->node_pt(n_node_1d - 1);
392
393 // The second node
394 node2_pt = quad_el_pt->node_pt((n_node_1d * n_node_1d) - 1);
395
396 // Break
397 break;
398 case QuadTreeNames::S:
399 // The first node
400 node1_pt = quad_el_pt->node_pt(0);
401
402 // The second node
403 node2_pt = quad_el_pt->node_pt(n_node_1d - 1);
404
405 // Break
406 break;
407 case QuadTreeNames::W:
408 // The first node
409 node1_pt = quad_el_pt->node_pt(n_node_1d * (n_node_1d - 1));
410
411 // The second node
412 node2_pt = quad_el_pt->node_pt(0);
413
414 // Break
415 break;
416 default:
417 // Create an ostringstream object to create an error message
418 std::ostringstream error_message_stream;
419
420 // Create the error message
421 error_message_stream
422 << "Input is " << i_edge << " but it can only\n"
423 << "be either N/E/S/W, i.e. " << QuadTreeNames::N << ","
424 << QuadTreeNames::E << "," << QuadTreeNames::S << " or "
425 << QuadTreeNames::W << std::endl;
426
427 // Throw an error
428 throw OomphLibError(error_message_stream.str(),
429 OOMPH_CURRENT_FUNCTION,
430 OOMPH_EXCEPTION_LOCATION);
431 }
432
433 // Create the edge associated with the i_edge-th edge of this
434 // element
435 Edge edge(node1_pt, node2_pt);
436
437 // Pair it up with the z-value at the current node slice
438 std::pair<Edge, double> edge_and_z_value_pair(edge, z_value);
439
440 // Store the edge (with its corresponding z-value)
441 edges_and_z_value.push_back(edge_and_z_value_pair);
442
443 // Is there is a matching entry in the map?
444 std::map<std::pair<Edge, double>, Vector<Node*>>::iterator it =
445 edge_to_nodes_map.find(edge_and_z_value_pair);
446
447 // If we're not at the end then it has already been created
448 if (it != edge_to_nodes_map.end())
449 {
450 // Get the vector of node pointers
451 Vector<Node*> edge_nodes_pt = it->second;
452
453#ifdef PARANOID
454 // Sanity check; make sure enough nodes have been stored
455 if (edge_nodes_pt.size() != n_node_1d)
456 {
457 // Throw an error
458 throw OomphLibError("Not enough nodes in edge_nodes_pt!",
459 OOMPH_CURRENT_FUNCTION,
460 OOMPH_EXCEPTION_LOCATION);
461 }
462#endif
463
464 // Remember that the i_edge-th edge has already been created!
465 switch (i_edge)
466 {
467 case QuadTreeNames::N:
468 // Indicate that the Northern edge has been done
469 has_edge_been_done[0] = true;
470
471 // Break
472 break;
473 case QuadTreeNames::E:
474 // Indicate that the Eastern edge has been done
475 has_edge_been_done[1] = true;
476
477 // Break
478 break;
479 case QuadTreeNames::S:
480 // Indicate that the Southern edge has been done
481 has_edge_been_done[2] = true;
482
483 // Break
484 break;
485 case QuadTreeNames::W:
486 // Indicate that the Western edge has been done
487 has_edge_been_done[3] = true;
488
489 // Break
490 break;
491 default:
492 // Throw an error
493 throw OomphLibError("Incorrect edge type!",
494 OOMPH_CURRENT_FUNCTION,
495 OOMPH_EXCEPTION_LOCATION);
496 }
497
498 //---------------------------------------------------------------
499 // Assign the nodes appropriately depending on which edge has
500 // already been created.
501 // NOTE: the nodes will have been created in anti-clockwise
502 // order which means we will have to assign them in clockwise
503 // order (from the perspective of the current element).
504 //---------------------------------------------------------------
505 // Loop over the nodes
506 for (unsigned i_node = 0; i_node < n_node_1d; i_node++)
507 {
508 // Calculate the value of quad_local_node_number depending on
509 // which edge the node lies on
510 switch (i_edge)
511 {
512 case QuadTreeNames::N:
513 // Calculate this Node's quad local node number
514 quad_local_node_number =
515 ((n_node_1d * n_node_1d) - 1) - i_node;
516
517 // Break
518 break;
519 case QuadTreeNames::E:
520 // Calculate this Node's quad local node number
521 quad_local_node_number =
522 (n_node_1d - 1) + (i_node * n_node_1d);
523
524 // Break
525 break;
526 case QuadTreeNames::S:
527 // Calculate this Node's quad local node number
528 quad_local_node_number = i_node;
529
530 // Break
531 break;
532 case QuadTreeNames::W:
533 // Calculate this Node's quad local node number
534 quad_local_node_number =
535 (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
536
537 // Break
538 break;
539 }
540
541 // Calculate this Node's cube local node number
542 cube_local_node_number =
543 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
544
545 // Assign the i_node-th node
546 cube_el_pt->node_pt(cube_local_node_number) =
547 edge_nodes_pt[(n_node_1d - 1) - i_node];
548
549 // Indicate that the appropriate node has been set up
550 has_node_been_setup[quad_local_node_number] = true;
551 } // for (unsigned i_node=0;i_node<n_node_1d;i_node++)
552 } // if (it!=edge_to_nodes_map.end())
553 } // for (unsigned i_edge=0;i<n_edge;i_edge++)
554
555 // The number of corners
556 unsigned n_corner = 4;
557
558 // Loop over the corner nodes
559 for (unsigned i_corner = 0; i_corner < n_corner; i_corner++)
560 {
561 // Calculate the value of quad_local_node_number depending on
562 // which corner the node lies on
563 switch (i_corner)
564 {
565 case 0:
566 // Bottom-left corner
567 quad_local_node_number = 0;
568
569 // Break
570 break;
571 case 1:
572 // Bottom-right corner
573 quad_local_node_number = n_node_1d - 1;
574
575 // Break
576 break;
577 case 2:
578 // Top-left corner
579 quad_local_node_number = n_node_1d * (n_node_1d - 1);
580
581 // Break
582 break;
583 case 3:
584 // Top-right corner
585 quad_local_node_number = (n_node_1d * n_node_1d) - 1;
586
587 // Break
588 break;
589 }
590
591 // If the node hasn't already been copied from some other edge
592 if (!has_node_been_setup[quad_local_node_number])
593 {
594 // Get a pointer to the first corner node (in the quad mesh)
595 Node* quad_corner_node_pt =
596 quad_el_pt->node_pt(quad_local_node_number);
597
598 // Pair up the corner node and z-value
599 std::pair<Node*, double> quad_corner_node_and_z_value(
600 quad_corner_node_pt, z_value);
601
602 // Have we made this Node yet?
603 std::map<std::pair<Node*, double>, Node*>::iterator it =
604 quad_corner_node_to_cube_node_map.find(
605 quad_corner_node_and_z_value);
606
607 // If we found a corresponding entry
608 if (it != quad_corner_node_to_cube_node_map.end())
609 {
610 // Index of the local node number
611 cube_local_node_number =
612 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
613
614 // Get the node and store it in the element
615 cube_el_pt->node_pt(cube_local_node_number) = it->second;
616
617 // Indicate that the node has now been created
618 has_node_been_setup[quad_local_node_number] = true;
619 }
620 } // if (!has_node_been_setup[quad_local_node_number])
621 } // for (unsigned i_corner=0;i_corner<n_corner;i_corner++)
622
623 //------------------------------------------------------------------
624 // By this point, all nodes which have already been created have
625 // been copied over. All the other nodes now need to be constructed.
626 //------------------------------------------------------------------
627 // Loop over the nodes in the second local coordinate direction
628 for (unsigned n_1 = 0; n_1 < n_node_1d; n_1++)
629 {
630 // Loop over the nodes in the first local coordinate direction
631 for (unsigned n_0 = 0; n_0 < n_node_1d; n_0++)
632 {
633 // Calculate the local node number (in the input mesh)
634 quad_local_node_number = n_0 + (n_node_1d * n_1);
635
636 // Calculate the local node number (in the extruded mesh)
637 cube_local_node_number =
638 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
639
640 // Check if the node has already been created (i.e. by copying)
641 if (has_node_been_setup[quad_local_node_number])
642 {
643 // If this node has already been set up then skip it
644 continue;
645 }
646
647 // Pointer to point to the node we are about to create
648 Node* cube_node_pt = 0;
649
650 //--------------------------------------------------------------
651 /// / Check if the node that we're about to construct lies on
652 /// the
653 // Front or back face of the mesh. The node can only lie on the
654 // front face if we're on the first element slice and first node
655 // slice. In contrast, the node can only lie on the back face if
656 // we're on the final element slice and final node slice.
657 //--------------------------------------------------------------
658 // The node lies on the front face or the back face
659 if (((i == 0) && (n_2 == 0)) ||
660 ((i == Nz - 1) && (n_2 == n_node_1d - 1)))
661 {
662 // The node lies on the front face
663 if ((i == 0) && (n_2 == 0))
664 {
665 // Create the boundary node
666 cube_node_pt = cube_el_pt->construct_boundary_node(
667 cube_local_node_number, time_stepper_pt);
668
669 // Indicate that the node has been created now
670 has_node_been_setup[quad_local_node_number] = true;
671
672 // Add the node to the appropriate boundary
673 add_boundary_node(n_boundary_in_quad_mesh, cube_node_pt);
674 }
675 // The node lies on the back face
676 else if ((i == Nz - 1) && (n_2 == n_node_1d - 1))
677 {
678 // Create the boundary node
679 cube_node_pt = cube_el_pt->construct_boundary_node(
680 cube_local_node_number, time_stepper_pt);
681
682 // Indicate that the node has been created now
683 has_node_been_setup[quad_local_node_number] = true;
684
685 // Add the node to the appropriate boundary
686 add_boundary_node(n_boundary_in_quad_mesh + 1,
687 cube_node_pt);
688 } // if ((i==0)&&(n_2==0))
689
690 // Get a pointer to the matching quad mesh node
691 Node* quad_node_pt =
692 quad_el_pt->node_pt(quad_local_node_number);
693
694 // If this Node lies on a boundary
695 if (quad_node_pt->is_on_boundary())
696 {
697 // Storage for the set of boundaries that this Node lies on
698 std::set<unsigned>* boundaries_pt;
699
700 // Get the boundaries that this Node lies on
701 quad_node_pt->get_boundaries_pt(boundaries_pt);
702
703 // Create an iterator to loop over the entries of the set
704 std::set<unsigned>::const_iterator const_it;
705
706 // Iterate over the boundaries
707 for (const_it = boundaries_pt->begin();
708 const_it != boundaries_pt->end();
709 const_it++)
710 {
711 // Add the node to the appropriate boundary (if it hasn't
712 // already been added)
713 add_boundary_node(*const_it, cube_node_pt);
714 }
715 } // quad_node_pt->is_on_boundary()
716 } // if (((i==0)&&(n_2==0))||((i==Nz-1)&&(n_2==n_node_1d-1)))
717
718
719 // Check if there is any other boundary information
720 if (el_boundary_information.size() > 0)
721 {
722 // How many boundaries are there to cover?
723 unsigned n_boundaries = el_boundary_information.size();
724
725 // Loop over the boundaries
726 for (unsigned i_bound = 0; i_bound < n_boundaries; i_bound++)
727 {
728 // If the node still hasn't been created yet
729 if (!has_node_been_setup[quad_local_node_number])
730 {
731 // Create the boundary node
732 cube_node_pt = cube_el_pt->construct_boundary_node(
733 cube_local_node_number, time_stepper_pt);
734
735 // Indicate that the node has been created now
736 has_node_been_setup[quad_local_node_number] = true;
737 }
738
739 // Get the boundary information associated with the
740 // i_bound-th entry
741 std::pair<unsigned, int> boundary_and_face_index =
742 el_boundary_information[i_bound];
743
744 // Decide whether or not the node lies on the boundary. If
745 // it does then add it to the boundary
746 switch (boundary_and_face_index.second)
747 {
748 case -1:
749 // If the s[0]=-1 boundary of the element lies on the
750 // boundary the node can only lie on the boundary if
751 // n_0=0.
752 if (n_0 == 0)
753 {
754 // Now add the node to the appropriate boundary
755 add_boundary_node(boundary_and_face_index.first,
756 cube_node_pt);
757 }
758
759 // Break
760 break;
761 case 1:
762 // If the s[0]=1 boundary of the element lies on the
763 // boundary the node can only lie on the boundary if
764 // n_0=n_node_1d-1
765 if (n_0 == n_node_1d - 1)
766 {
767 // Now add the node to the appropriate boundary
768 add_boundary_node(boundary_and_face_index.first,
769 cube_node_pt);
770 }
771
772 // Break
773 break;
774 case -2:
775 // If the s[1]=-1 boundary of the element lies on the
776 // boundary the node can only lie on the boundary if
777 // n_1=0.
778 if (n_1 == 0)
779 {
780 // Now add the node to the appropriate boundary
781 add_boundary_node(boundary_and_face_index.first,
782 cube_node_pt);
783 }
784
785 // Break
786 break;
787 case 2:
788 // If the s[1]=1 boundary of the element lies on the
789 // boundary the node can only lie on the boundary if
790 // n_1=n_node_1d-1.
791 if (n_1 == n_node_1d - 1)
792 {
793 // Now add the node to the appropriate boundary
794 add_boundary_node(boundary_and_face_index.first,
795 cube_node_pt);
796 }
797
798 // Break
799 break;
800 } // switch (boundary_and_face_index.second)
801 } // for (unsigned i_bound=0;i_bound<n_boundaries;i_bound++)
802 } // if (el_boundary_information.size()>0)
803
804 // Okay, so this Node does not lie on a face of a quad element
805 // that lies on a boundary. However(!), there is still another
806 // case in which it might lie on a mesh boundary; when the Node
807 // is the ONLY Node in the element that lies on the boundary.
808 // It's a slightly odd scenario but can easily occur on meshes
809 // created using the QuadFromTriangleMesh class.
810 Node* quad_node_pt =
811 quad_el_pt->node_pt(quad_local_node_number);
812
813 // If this Node lies on a boundary
814 if (quad_node_pt->is_on_boundary())
815 {
816 // If we haven't created the Node yet AND the corresponding
817 // Node in the quad mesh definitely lies on the boundary
818 if (cube_node_pt == 0)
819 {
820 // Create the boundary node
821 cube_node_pt = cube_el_pt->construct_boundary_node(
822 cube_local_node_number, time_stepper_pt);
823
824 // Indicate that the node has been created now
825 has_node_been_setup[quad_local_node_number] = true;
826 }
827
828 // Storage for the set of boundaries that this Node lies on
829 std::set<unsigned>* boundaries_pt;
830
831 // Get the boundaries that this Node lies on
832 quad_node_pt->get_boundaries_pt(boundaries_pt);
833
834 // Create an iterator to loop over the entries of the set
835 std::set<unsigned>::const_iterator const_it;
836
837 // Iterate over the boundaries
838 for (const_it = boundaries_pt->begin();
839 const_it != boundaries_pt->end();
840 const_it++)
841 {
842 // Add the node to the appropriate boundary (if it hasn't
843 // already been added)
844 add_boundary_node(*const_it, cube_node_pt);
845 }
846 } // quad_node_pt->is_on_boundary()
847
848 // If the node still hasn't been created yet
849 if (!has_node_been_setup[quad_local_node_number])
850 {
851 // The only other possibility is that it is a regular node
852 cube_node_pt = cube_el_pt->construct_node(
853 cube_local_node_number, time_stepper_pt);
854
855 // Indicate that the node has now been set up
856 has_node_been_setup[quad_local_node_number] = true;
857 }
858
859 // Set the x-coordinate of the node
860 cube_node_pt->x(0) = quad_node_pt->x(0);
861
862 // Set the y-coordinate of the node
863 cube_node_pt->x(1) = quad_node_pt->x(1);
864
865 // Set the z-coordinate of the node
866 cube_node_pt->x(2) = z_value;
867
868 // Add it to global storage
869 Node_pt[node_count] = cube_node_pt;
870
871 // Increment the counter
872 node_count++;
873 } // for (unsigned n_0=0;n_0<n_node_1d;n_0++)
874 } // for (unsigned n_1=0;n_1<n_node_1d;n_1++)
875
876 //------------------------------------------------------------------
877 // We've set up the nodes in one slice, now set up the edge-to-node
878 // information and store it in the map edge_to_nodes_map.
879 //------------------------------------------------------------------
880 for (unsigned i_edge = 0; i_edge < n_edge; i_edge++)
881 {
882 // If it hasn't already been done
883 if (!has_edge_been_done[i_edge])
884 {
885 // Create a vector to store the nodes on this edge
886 Vector<Node*> edge_nodes_pt(n_node_1d, 0);
887
888 // Loop over the nodes and add them
889 for (unsigned i_node = 0; i_node < n_node_1d; i_node++)
890 {
891 // Calculate the value of quad_local_node_number depending on
892 // which corner the node lies on
893 switch (i_edge)
894 {
895 case 0:
896 // Northern edge
897 quad_local_node_number =
898 ((n_node_1d * n_node_1d) - 1) - i_node;
899
900 // Break
901 break;
902 case 1:
903 // Eastern edge
904 quad_local_node_number =
905 (n_node_1d - 1) + (i_node * n_node_1d);
906
907 // Break
908 break;
909 case 2:
910 // Southern edge
911 quad_local_node_number = i_node;
912
913 // Break
914 break;
915 case 3:
916 // Western edge
917 quad_local_node_number =
918 (n_node_1d * (n_node_1d - 1)) - (i_node * n_node_1d);
919
920 // Break
921 break;
922 }
923
924 // Index of the local node number
925 cube_local_node_number =
926 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
927
928 // Store the i_node-th node along this edge
929 edge_nodes_pt[i_node] =
930 cube_el_pt->node_pt(cube_local_node_number);
931 } // for (unsigned i_node=0;i_node<n_node_1d;i_node++)
932
933 // Store the nodes along the north edge of the element
934 edge_to_nodes_map[edges_and_z_value[i_edge]] = edge_nodes_pt;
935 } // if (!has_edge_been_done[i_edge])
936 } // for (unsigned i_edge=0;i_edge<n_edge;i_edge++)
937
938 //--------------------------------------------------------------------
939 // Now store the corner nodes of the element at the current node
940 // slice in the map quad_corner_node_to_cube_node_map.
941 //--------------------------------------------------------------------
942 // Loop over the corner nodes
943 for (unsigned i_corner = 0; i_corner < n_corner; i_corner++)
944 {
945 // Calculate the value of quad_local_node_number depending on
946 // which corner the node lies on
947 switch (i_corner)
948 {
949 case 0:
950 // Bottom-left corner
951 quad_local_node_number = 0;
952
953 // Break
954 break;
955 case 1:
956 // Bottom-right corner
957 quad_local_node_number = n_node_1d - 1;
958
959 // Break
960 break;
961 case 2:
962 // Top-left corner
963 quad_local_node_number = n_node_1d * (n_node_1d - 1);
964
965 // Break
966 break;
967 case 3:
968 // Top-right corner
969 quad_local_node_number = (n_node_1d * n_node_1d) - 1;
970
971 // Break
972 break;
973 }
974
975 // Index of the corresponding corner node (in the cube element)
976 cube_local_node_number =
977 quad_local_node_number + (n_node_1d * n_node_1d * n_2);
978
979 // Get a pointer to the first corner node (in the quad mesh)
980 Node* quad_corner_node_pt =
981 quad_el_pt->node_pt(quad_local_node_number);
982
983 // Get a pointer to the first corner node (in the cube mesh)
984 Node* cube_corner_node_pt =
985 cube_el_pt->node_pt(cube_local_node_number);
986
987 // Pair up the corner node and z-value
988 std::pair<Node*, double> quad_corner_node_and_z_value(
989 quad_corner_node_pt, z_value);
990
991 // Store it in the map
992 quad_corner_node_to_cube_node_map[quad_corner_node_and_z_value] =
993 cube_corner_node_pt;
994 } // for (unsigned i_corner=0;i_corner<n_corner;i_corner++)
995
996#ifdef PARANOID
997 // Loop over all the nodes in the current node slice
998 for (unsigned i_node = 0; i_node < (n_node_1d * n_node_1d);
999 i_node++)
1000 {
1001 // Make sure every node in the current slice has been set up
1002 if (!has_node_been_setup[i_node])
1003 {
1004 // Create an ostringstream object to create an error message
1005 std::ostringstream error_message_stream;
1006
1007 // Create the error message
1008 error_message_stream << "There are nodes in element " << e
1009 << " which have not been constructed!"
1010 << std::endl;
1011
1012 // Throw an error
1013 throw OomphLibError(error_message_stream.str(),
1014 OOMPH_CURRENT_FUNCTION,
1015 OOMPH_EXCEPTION_LOCATION);
1016 }
1017 } // for (unsigned i_node=0;i_node<(n_node_1d*n_node_1d);i_node++)
1018#endif
1019 } // if ((i>0)&&(n_2==0))
1020 } // for (unsigned n_2=0;n_2<n_node_1d;n_2++)
1021 } // for (unsigned e=0;e<n_element_in_quad_mesh;e++)
1022
1023 // Sanity check; make sure enough nodes have been stored
1024 if (node_count != (1 + (n_node_1d - 1) * (i + 1)) * n_node_in_quad_mesh)
1025 {
1026 // Create an ostringstream object to create an error message
1027 std::ostringstream error_message_stream;
1028
1029 // The number of nodes we expect there to be in the mesh
1030 int expected_n_node =
1031 (1 + (n_node_1d - 1) * (i + 1)) * n_node_in_quad_mesh;
1032
1033 // Difference in node count
1034 int node_diff = expected_n_node - node_count;
1035
1036 // Create the error message
1037 error_message_stream
1038 << "There are meant to be " << expected_n_node
1039 << " nodes in the extruded mesh by the\nend of slice " << i
1040 << " but you have " << node_count << ". "
1041 << "That's a difference of " << node_diff << "!";
1042
1043 // Throw an error
1044 throw OomphLibError(error_message_stream.str(),
1045 OOMPH_CURRENT_FUNCTION,
1046 OOMPH_EXCEPTION_LOCATION);
1047 }
1048 } // for (unsigned i=0;i<Nz;i++)
1049
1050 //------------------------------------------------------------------------
1051 // Count the number of nodes on the boundaries of the quad mesh; there
1052 // should really be a helper function that does this (hint hint)
1053 //------------------------------------------------------------------------
1054 // Use a set to make sure we don't count edge/corner nodes several times
1055 std::set<const Node*> all_quad_boundary_nodes_pt;
1056
1057 // The number of boundary nodes in our extruded mesh
1058 int n_quad_mesh_boundary_node = 0;
1059
1060 // The number of boundaries in the quad mesh
1061 for (unsigned b = 0; b < n_boundary_in_quad_mesh; b++)
1062 {
1063 // Add on the number of nodes on the b-th boundary
1064 unsigned n_boundary_node = quad_mesh_pt->nboundary_node(b);
1065
1066 // Loop over the nodes on this boundary
1067 for (unsigned n = 0; n < n_boundary_node; n++)
1068 {
1069 // We haven't come across this boundary node yet
1070 if (all_quad_boundary_nodes_pt.find(quad_mesh_pt->boundary_node_pt(
1071 b, n)) == all_quad_boundary_nodes_pt.end())
1072 {
1073 // Count this node
1074 n_quad_mesh_boundary_node++;
1075
1076 // Now remember it so we don't count it twice
1077 all_quad_boundary_nodes_pt.insert(
1078 quad_mesh_pt->boundary_node_pt(b, n));
1079 }
1080 } // for (unsigned n=0; n<n_boundary_node; n++)
1081 } // for (unsigned b=0; b<n_boundary_in_quad_mesh; b++)
1082
1083 //------------------------------------------------------------------------
1084 // Count the number of nodes on the boundaries of the extruded cube mesh;
1085 // (recall hint hint from above...)
1086 //------------------------------------------------------------------------
1087 // Number of boundary nodes expected in the extruded mesh; initialise to
1088 // the number of boundary nodes on the front and back faces
1089 int n_expected_boundary_node = 2 * quad_mesh_pt->nnode();
1090
1091 // Add on the nodes expected on the extruded mesh boundaries (but ignoring
1092 // the boundary nodes on the front and back faces)
1093 n_expected_boundary_node +=
1094 (((n_node_1d - 1) * Nz) - 1) * n_quad_mesh_boundary_node;
1095
1096 // The number of boundary nodes in our extruded mesh
1097 int n_extruded_mesh_boundary_node = 0;
1098
1099 // Use a set to make sure we don't count edge/corner nodes several times
1100 std::set<const Node*> all_extruded_boundary_nodes_pt;
1101
1102 // Loop over the boundaries of the mesh
1103 for (unsigned b = 0; b < n_boundary_in_extruded_mesh; b++)
1104 {
1105 // Add on the number of nodes on the b-th boundary
1106 unsigned n_boundary_node = nboundary_node(b);
1107
1108 // Loop over the nodes on this boundary
1109 for (unsigned n = 0; n < n_boundary_node; n++)
1110 {
1111 // We haven't come across this boundary node yet
1112 if (all_extruded_boundary_nodes_pt.find(this->boundary_node_pt(b, n)) ==
1113 all_extruded_boundary_nodes_pt.end())
1114 {
1115 // Count this node
1116 n_extruded_mesh_boundary_node++;
1117
1118 // Now remember it so we don't count it twice
1119 all_extruded_boundary_nodes_pt.insert(this->boundary_node_pt(b, n));
1120 }
1121 } // for (unsigned n=0; n<n_boundary_node; n++)
1122 } // for (unsigned b=0; b<n_boundary_in_extruded_mesh; b++)
1123
1124 //------------------------------------------------------------------------
1125 // Check that there is the correct number of boundary nodes
1126 //------------------------------------------------------------------------
1127 // Error: we don't have the correct number of boundary nodes
1128 if (n_extruded_mesh_boundary_node != n_expected_boundary_node)
1129 {
1130 // Create an ostringstream object to create an error message
1131 std::ostringstream error_message_stream;
1132
1133 // Difference in node count
1134 int node_diff = n_expected_boundary_node - n_extruded_mesh_boundary_node;
1135
1136 // Create the error message
1137 error_message_stream << "There should be " << n_expected_boundary_node
1138 << " boundary nodes in the extruded mesh but there"
1139 << "\nare only " << n_extruded_mesh_boundary_node
1140 << " boundary nodes. That's a difference of "
1141 << node_diff << "!" << std::endl;
1142
1143 // Throw an error
1144 throw OomphLibError(error_message_stream.str(),
1145 OOMPH_CURRENT_FUNCTION,
1146 OOMPH_EXCEPTION_LOCATION);
1147 }
1148
1149 // If the user wishes the mesh setup time to be doc-ed
1150 if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1151 {
1152 // Tell the user
1153 oomph_info << "Time taken to extrude mesh [sec]: "
1154 << TimingHelpers::timer() - t_mesh_construct_start
1155 << std::endl;
1156 }
1157
1158 // Get the current (i.e. start) time
1159 double t_mesh_macro_start = TimingHelpers::timer();
1160
1161 //----------------------------------------------------------------
1162 // If an element has a macro-element representation in the quad
1163 // mesh then it need an extruded macro-element representation in
1164 // the extruded mesh so set that up here. Then, use the extruded
1165 // macro-element representation of the elements to move all the
1166 // nodes that need to be moved.
1167 //----------------------------------------------------------------
1168 // Create a map that takes a Domain and returns the corresponding
1169 // ExtrudedDomain object
1170 std::map<Domain*, ExtrudedDomain*> domain_to_extruded_domain_map;
1171
1172 // Loop over the elements in the input mesh
1173 for (unsigned e = 0; e < n_element_in_quad_mesh; e++)
1174 {
1175 // Get a pointer to the e-th element in the input mesh.
1176 // NOTE: The element is upcast to the QElementBase class so we
1177 // can access the s_macro_ll() and s_macro_ur() functions.
1178 QElementBase* quad_el_pt =
1179 dynamic_cast<QElementBase*>(quad_mesh_pt->finite_element_pt(e));
1180
1181 // Check if the element has a macro-element representation
1182 if (quad_el_pt->macro_elem_pt() != 0)
1183 {
1184 // Get a pointer to the corresponding Domain object
1185 Domain* domain_pt = quad_el_pt->macro_elem_pt()->domain_pt();
1186
1187 // How many macro elements are there in this Domain?
1188 unsigned n_macro_element = domain_pt->nmacro_element();
1189
1190 // Create a pointer to point to the extruded version of the above Domain
1191 ExtrudedDomain* extruded_domain_pt = 0;
1192
1193 // Create an iterator for the map
1194 std::map<Domain*, ExtrudedDomain*>::iterator it;
1195
1196 // Check if there is a corresponding entry in the map
1197 it = domain_to_extruded_domain_map.find(domain_pt);
1198
1199 // Check if the extruded version of this domain has already been created
1200 if (it != domain_to_extruded_domain_map.end())
1201 {
1202 // Get the associated ExtrudedDomain pointer
1203 extruded_domain_pt = it->second;
1204 }
1205 // If there is no entry then we need to create the extruded domain
1206 else
1207 {
1208 // Create an ExtrudedDomain object associated with the mesh
1209 extruded_domain_pt = new ExtrudedDomain(domain_pt, Nz, Zmin, Zmax);
1210
1211 // Now store it in the map
1212 domain_to_extruded_domain_map[domain_pt] = extruded_domain_pt;
1213
1214 // Add it to the mesh's list of extruded domains
1215 Extruded_domain_pt.push_back(extruded_domain_pt);
1216 }
1217
1218 // Loop over each element slice in the z-direction
1219 for (unsigned i = 0; i < Nz; i++)
1220 {
1221 // What is the ID number of the macro element associated with
1222 // this quad element?
1223 unsigned macro_id =
1224 (i * n_macro_element +
1225 quad_el_pt->macro_elem_pt()->macro_element_number());
1226
1227 // Get a pointer to the extruded macro element
1228 ExtrudedMacroElement* extruded_macro_element_pt =
1229 extruded_domain_pt->macro_element_pt(macro_id);
1230
1231 // Get the index of the cube element in the extruded mesh
1232 unsigned element_index = i * n_element_in_quad_mesh + e;
1233
1234 // Get the corresponding extruded element from the global storage
1235 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Element_pt[element_index]);
1236
1237 // Set its macro element pointer
1238 el_pt->set_macro_elem_pt(extruded_macro_element_pt);
1239
1240 // The number of dimensions in the extruded mesh
1241 unsigned n_dim = 3;
1242
1243 // Loop over the coordinate directions
1244 for (unsigned i = 0; i < n_dim - 1; i++)
1245 {
1246 // The i-th local coordinate value of the "lower-left" corner of
1247 // the current element
1248 el_pt->s_macro_ll(i) = quad_el_pt->s_macro_ll(i);
1249
1250 // The i-th local coordinate value of the "upper-right" corner of
1251 // the current element
1252 el_pt->s_macro_ur(i) = quad_el_pt->s_macro_ur(i);
1253 }
1254
1255 // The local coordinate value of the "lower-left" corner which
1256 // corresponds to the time direction
1257 el_pt->s_macro_ll(n_dim - 1) = -1.0;
1258
1259 // The local coordinate value of the "upper-right" corner which
1260 // corresponds to the time direction
1261 el_pt->s_macro_ur(n_dim - 1) = 1.0;
1262 } // for (unsigned i=0;i<Nz;i++)
1263 } // if (quad_el_pt->macro_element_pt()!=0)
1264 } // for (unsigned e=0;e<n_element_in_quad_mesh;e++)
1265
1266 // If the user wishes the mesh setup time to be doc-ed
1267 if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1268 {
1269 // Tell the user
1270 oomph_info << "Time taken to set up extruded macro-elements [sec]: "
1271 << TimingHelpers::timer() - t_mesh_macro_start << std::endl;
1272 }
1273
1274 // Call the node update function to do all the relevant node moving
1275 node_update();
1276
1277#ifdef PARANOID
1278 // Make sure there are no duplicate nodes (i.e. two or more nodes that
1279 // occupy the same Eulerian position)
1280 if (check_for_repeated_nodes())
1281 {
1282 // Throw a warning
1283 OomphLibWarning("The mesh contains repeated nodes!",
1284 OOMPH_CURRENT_FUNCTION,
1285 OOMPH_EXCEPTION_LOCATION);
1286 }
1287
1288 // Also make sure enough nodes have been been created
1289 if (node_count != n_node_in_cube_mesh)
1290 {
1291 // Create an ostringstream object to create an error message
1292 std::ostringstream error_message_stream;
1293
1294 // Create the error message
1295 error_message_stream << "An incorrect number of nodes have been created! "
1296 << "There are\nmeant to be " << n_node_in_cube_mesh
1297 << " nodes in the extruded mesh but\nonly "
1298 << node_count << " have actually been constructed!"
1299 << std::endl;
1300
1301 // Throw a warning
1302 OomphLibWarning(error_message_stream.str(),
1303 OOMPH_CURRENT_FUNCTION,
1304 OOMPH_EXCEPTION_LOCATION);
1305 }
1306#endif
1307
1308 // Setup lookup scheme that establishes which elements are located
1309 // on the mesh boundaries
1310 setup_boundary_element_info();
1311
1312 // If the user wishes the mesh setup time to be doc-ed
1313 if (MeshExtrusionHelpers::Mesh_extrusion_helper.doc_mesh_setup_time())
1314 {
1315 // Tell the user
1316 oomph_info << "Total time for extruded mesh creation/setup [sec]: "
1317 << TimingHelpers::timer() - t_mesh_setup_start << std::endl;
1318 }
1319 } // End of build_mesh
1320} // End of namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:67
unsigned nmacro_element()
Number of macro elements in domain.
Definition: domain.h:123
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
Definition: mesh.h:2692
void build_mesh(QuadMeshBase *quad_mesh_pt, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Generic mesh construction function: contains all the hard work.
Base class for ExtrudedDomains with curvilinear and/or time-dependent boundaries. ExtrudedDomain boun...
ExtrudedMacroElement * macro_element_pt(const unsigned &i)
Access to i-th extruded macro element.
DRAIG: FILL IN COMPLETE DESCRIPTION ONCE FINISHED...
A general Finite Element class.
Definition: elements.h:1313
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:2218
Domain *& domain_pt()
Access function to the Domain_pt.
unsigned & macro_element_number()
Access function to the Macro_element_number.
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
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:1365
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition: Qelements.h:91
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:202
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:186
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: quad_mesh.h:57
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
class oomph::MeshExtrusionHelpers::ExtrusionHelper Mesh_extrusion_helper
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...