simple_cubic_scaffold_tet_mesh.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2025 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#include "Telements.h"
28
29
30namespace oomph
31{
32 //===========================================================
33 /// Scaffold mesh for cubic tet mesh.
34 //===========================================================
36 const unsigned& n_x,
37 const unsigned& n_y,
38 const unsigned& n_z,
39 const double& l_x,
40 const double& l_y,
41 const double& l_z,
42 TimeStepper* time_stepper_pt)
43 {
44 // Storage for the pointers to the nodes
45 RankThreeTensor<Node*> vertex_node_pt(n_x + 1, n_y + 1, n_z + 1);
50
51 // "Lower left" corner
52 double x_0 = 0.0;
53 double y_0 = 0.0;
54 double z_0 = 0.0;
55
56 // Increments
57 double dx = l_x / double(n_x);
58 double dy = l_y / double(n_y);
59 double dz = l_z / double(n_z);
60
61 // Total number of nodes
62 unsigned nnod = (n_x + 1) * (n_y + 1) * (n_z + 1) + // vertex nodes
63 n_x * n_y * n_z + // central nodes
64 (n_x * n_y) * (n_z + 1) + // faces
65 (n_x * n_z) * (n_y + 1) + (n_y * n_z) * (n_x + 1);
66 Node_pt.reserve(nnod);
67
68 // Total number of elements
69 unsigned nelem = 24 * n_x * n_y * n_z;
70 Element_pt.reserve(nelem);
71
72 // Six boundaries
74
75
76 // Generate the vertex nodes of all cells
77 //=======================================
79 for (unsigned k = 0; k < n_z + 1; k++)
80 {
81 for (unsigned j = 0; j < n_y + 1; j++)
82 {
83 for (unsigned i = 0; i < n_x + 1; i++)
84 {
85 // Need to work out whether to create a boundary node or not
86 if (((i < n_x + 1) && (j < n_y + 1) && (k < n_z + 1)) &&
87 ((i == 0) || (i == n_x) || (j == 0) || (j == n_y) || (k == 0) ||
88 (k == n_z)))
89 {
90 node_pt = new BoundaryNode<Node>(3, 1, 0);
91 }
92 else
93 {
94 node_pt = new Node(3, 1, 0);
95 }
96
97 vertex_node_pt(i, j, k) = node_pt;
98
99 // Coordinates
100 node_pt->x(0) = x_0 + double(i) * dx;
101 node_pt->x(1) = y_0 + double(j) * dy;
102 node_pt->x(2) = z_0 + double(k) * dz;
103
104 // Real node?
105 if ((i < n_x + 1) && (j < n_y + 1) && (k < n_z + 1))
106 {
107 // Add node
108 Node_pt.push_back(node_pt);
109
110 // Left boundary is 0
111 if (i == 0)
112 {
114 }
115 // Right boundary is 1
116 else if (i == n_x)
117 {
119 }
120
121 // Front boundary is 2
122 if (j == 0)
123 {
125 }
126 // Back boundary is 3
127 else if (j == n_y)
128 {
130 }
131
132 // Down boundary is 4
133 if (k == 0)
134 {
136 }
137 // Up boundary is 5
138 else if (k == n_z)
139 {
141 }
142 }
143 }
144 }
145 }
146
147
148 // Generate the face nodes of all cells
149 //=====================================
150 for (unsigned k = 0; k < n_z + 1; k++)
151 {
152 for (unsigned j = 0; j < n_y + 1; j++)
153 {
154 for (unsigned i = 0; i < n_x + 1; i++)
155 {
156 // Node on front face of cell
157 //---------------------------
158 // Need to work out whether to create boundary node or not
159 if (((i < n_x) && (k < n_z)) && ((j == 0) || (j == n_y)))
160 {
161 node_pt = new BoundaryNode<Node>(3, 1, 0);
162 }
163 else
164 {
165 node_pt = new Node(3, 1, 0);
166 }
168
169
170 // Coordinates
171 node_pt->x(0) = x_0 + 0.5 * dx + double(i) * dx;
172 node_pt->x(1) = y_0 + double(j) * dy;
173 node_pt->x(2) = z_0 + 0.5 * dz + double(k) * dz;
174
175 // Real node?
176 if ((i < n_x) && (k < n_z))
177 {
178 // Add node to mesh
179 Node_pt.push_back(node_pt);
180
181 // Front boundary is 2
182 if (j == 0)
183 {
185 }
186 // Back boundary is 3
187 else if (j == n_y)
188 {
190 }
191 }
192
193 // Node on left face of cell
194 //--------------------------
195 // Work out whether to construct boundary node or not
196 if (((j < n_y) && (k < n_z)) && ((i == 0) || (i == n_x)))
197 {
198 node_pt = new BoundaryNode<Node>(3, 1, 0);
199 }
200 else
201 {
202 node_pt = new Node(3, 1, 0);
203 }
205
206 // Coordinates
207 node_pt->x(0) = x_0 + double(i) * dx;
208 node_pt->x(1) = y_0 + 0.5 * dy + double(j) * dy;
209 node_pt->x(2) = z_0 + 0.5 * dz + double(k) * dz;
210
211
212 // Real node?
213 if ((j < n_y) && (k < n_z))
214 {
215 // Add node to mesh
216 Node_pt.push_back(node_pt);
217
218 // Left boundary is 0
219 if (i == 0)
220 {
222 }
223 // Right boundary is 1
224 else if (i == n_x)
225 {
227 }
228 }
229
230 // Node on down face of cell
231 //--------------------------
232 if (((i < n_x) && (j < n_y)) && ((k == 0) || (k == n_z)))
233 {
234 node_pt = new BoundaryNode<Node>(3, 1, 0);
235 }
236 else
237 {
238 node_pt = new Node(3, 1, 0);
239 }
241
242 // Coordinates
243 node_pt->x(0) = x_0 + 0.5 * dx + double(i) * dx;
244 node_pt->x(1) = y_0 + 0.5 * dy + double(j) * dy;
245 node_pt->x(2) = z_0 + double(k) * dz;
246
247
248 // Real node?
249 if ((i < n_x) && (j < n_y))
250 {
251 // Add node to mesh
252 Node_pt.push_back(node_pt);
253
254 // Down boundary is 4
255 if (k == 0)
256 {
258 }
259 // Up boundary is 5
260 else if (k == n_z)
261 {
263 }
264 }
265 }
266 }
267 }
268
269
270 // Central nodes for all cells
271 for (unsigned k = 0; k < n_z; k++)
272 {
273 for (unsigned j = 0; j < n_y; j++)
274 {
275 for (unsigned i = 0; i < n_x; i++)
276 {
277 node_pt = new Node(3, 1, 0);
279 Node_pt.push_back(node_pt);
280 node_pt->x(0) = x_0 + 0.5 * dx + double(i) * dx;
281 node_pt->x(1) = y_0 + 0.5 * dy + double(j) * dy;
282 node_pt->x(2) = z_0 + 0.5 * dz + double(k) * dz;
283 }
284 }
285 }
286
287
288 // Loop over blocks and create elements
290 for (unsigned k = 0; k < n_z; k++)
291 {
292 for (unsigned j = 0; j < n_y; j++)
293 {
294 for (unsigned i = 0; i < n_x; i++)
295 {
296 // FRONT FACE
297 //===========
298
299 // Left element on front face
300 //---------------------------
301 el_pt = new TElement<3, 2>;
302
303 // LFD
304 el_pt->node_pt(0) = vertex_node_pt(i, j, k);
305
306 // Central front face
308
309 // LFU
310 el_pt->node_pt(2) = vertex_node_pt(i, j, k + 1);
311
312 // central node
313 el_pt->node_pt(3) = central_node_pt(i, j, k);
314
315 Element_pt.push_back(el_pt);
316
317
318 // Down element on front face
319 //---------------------------
320 el_pt = new TElement<3, 2>;
321
322 // LFD
323 el_pt->node_pt(0) = vertex_node_pt(i, j, k);
324
325 // RFD
326 el_pt->node_pt(1) = vertex_node_pt(i + 1, j, k);
327
328 // Central front face
330
331 // Central node
332 el_pt->node_pt(3) = central_node_pt(i, j, k);
333
334 Element_pt.push_back(el_pt);
335
336
337 // Right element on front face
338 //---------------------------
339 el_pt = new TElement<3, 2>;
340
341 // RFD
342 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k);
343
344 // RFU
345 el_pt->node_pt(1) = vertex_node_pt(i + 1, j, k + 1);
346
347 // Central front face
349
350 // Central node
351 el_pt->node_pt(3) = central_node_pt(i, j, k);
352
353 Element_pt.push_back(el_pt);
354
355
356 // Up element on front face
357 //---------------------------
358 el_pt = new TElement<3, 2>;
359
360 // RFU
361 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k + 1);
362
363 // LFU
364 el_pt->node_pt(1) = vertex_node_pt(i, j, k + 1);
365
366 // Central front face
368
369 // Central node
370 el_pt->node_pt(3) = central_node_pt(i, j, k);
371
372 Element_pt.push_back(el_pt);
373
374
375 // RIGHT FACE
376 //===========
377
378 // Front element on right face
379 //---------------------------
380 el_pt = new TElement<3, 2>;
381
382 // RFD
383 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k);
384
385 // Central right face
386 el_pt->node_pt(1) = left_face_node_pt(i + 1, j, k);
387
388 // RFU
389 el_pt->node_pt(2) = vertex_node_pt(i + 1, j, k + 1);
390
391 // Central node
392 el_pt->node_pt(3) = central_node_pt(i, j, k);
393
394 Element_pt.push_back(el_pt);
395
396
397 // Down element on right face
398 //---------------------------
399 el_pt = new TElement<3, 2>;
400
401 // RFD
402 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k);
403
404 // RBD
405 el_pt->node_pt(1) = vertex_node_pt(i + 1, j + 1, k);
406
407 // Central right face
408 el_pt->node_pt(2) = left_face_node_pt(i + 1, j, k);
409
410 // central node
411 el_pt->node_pt(3) = central_node_pt(i, j, k);
412
413 Element_pt.push_back(el_pt);
414
415
416 // Back element on right face
417 //---------------------------
418 el_pt = new TElement<3, 2>;
419
420 // RBD
421 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k);
422
423 // RBU
424 el_pt->node_pt(1) = vertex_node_pt(i + 1, j + 1, k + 1);
425
426 // Central right face
427 el_pt->node_pt(2) = left_face_node_pt(i + 1, j, k);
428
429 // central node
430 el_pt->node_pt(3) = central_node_pt(i, j, k);
431
432 Element_pt.push_back(el_pt);
433
434
435 // Up element on right face
436 //---------------------------
437 el_pt = new TElement<3, 2>;
438
439 // RBU
440 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k + 1);
441
442 // RFU
443 el_pt->node_pt(1) = vertex_node_pt(i + 1, j, k + 1);
444
445 // Central right face
446 el_pt->node_pt(2) = left_face_node_pt(i + 1, j, k);
447
448 // Central node
449 el_pt->node_pt(3) = central_node_pt(i, j, k);
450
451 Element_pt.push_back(el_pt);
452
453
454 // UP FACE
455 //===========
456
457 // Front element on up face
458 //---------------------------
459 el_pt = new TElement<3, 2>;
460
461 // LFU
462 el_pt->node_pt(0) = vertex_node_pt(i, j, k + 1);
463
464 // RFU
465 el_pt->node_pt(1) = vertex_node_pt(i + 1, j, k + 1);
466
467 // Central up face
468 el_pt->node_pt(2) = down_face_node_pt(i, j, k + 1);
469
470 // Central node
471 el_pt->node_pt(3) = central_node_pt(i, j, k);
472
473 Element_pt.push_back(el_pt);
474
475
476 // Front element on up face
477 //---------------------------
478 el_pt = new TElement<3, 2>;
479
480 // RFU
481 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k + 1);
482
483 // RBU
484 el_pt->node_pt(1) = vertex_node_pt(i + 1, j + 1, k + 1);
485
486 // Central up face
487 el_pt->node_pt(2) = down_face_node_pt(i, j, k + 1);
488
489 // central node
490 el_pt->node_pt(3) = central_node_pt(i, j, k);
491
492 Element_pt.push_back(el_pt);
493
494
495 // Back element on up face
496 //---------------------------
497 el_pt = new TElement<3, 2>;
498
499 // RBU
500 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k + 1);
501
502 // LBU
503 el_pt->node_pt(1) = vertex_node_pt(i, j + 1, k + 1);
504
505 // Central up face
506 el_pt->node_pt(2) = down_face_node_pt(i, j, k + 1);
507
508 // central node
509 el_pt->node_pt(3) = central_node_pt(i, j, k);
510
511 Element_pt.push_back(el_pt);
512
513
514 // Left element on up face
515 //---------------------------
516 el_pt = new TElement<3, 2>;
517
518
519 // LBU
520 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k + 1);
521
522 // RFU
523 el_pt->node_pt(1) = vertex_node_pt(i, j, k + 1);
524
525 // Central up face
526 el_pt->node_pt(2) = down_face_node_pt(i, j, k + 1);
527
528 // Central node
529 el_pt->node_pt(3) = central_node_pt(i, j, k);
530
531 Element_pt.push_back(el_pt);
532
533
534 // DOWN FACE
535 //===========
536
537 // Front element on down face
538 //---------------------------
539 el_pt = new TElement<3, 2>;
540
541 // LFD
542 el_pt->node_pt(0) = vertex_node_pt(i, j, k);
543
544 // RFD
545 el_pt->node_pt(2) = vertex_node_pt(i + 1, j, k);
546
547 // Central down face
549
550 // Central node
551 el_pt->node_pt(3) = central_node_pt(i, j, k);
552
553 Element_pt.push_back(el_pt);
554
555
556 // Front element on down face
557 //---------------------------
558 el_pt = new TElement<3, 2>;
559
560 // RFD
561 el_pt->node_pt(0) = vertex_node_pt(i + 1, j, k);
562
563 // RBD
564 el_pt->node_pt(2) = vertex_node_pt(i + 1, j + 1, k);
565
566 // Central down face
568
569 // central node
570 el_pt->node_pt(3) = central_node_pt(i, j, k);
571
572 Element_pt.push_back(el_pt);
573
574
575 // Back element on down face
576 //---------------------------
577 el_pt = new TElement<3, 2>;
578
579 // RBD
580 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k);
581
582 // LBD
583 el_pt->node_pt(2) = vertex_node_pt(i, j + 1, k);
584
585 // Central down face
587
588 // central node
589 el_pt->node_pt(3) = central_node_pt(i, j, k);
590
591 Element_pt.push_back(el_pt);
592
593
594 // Left element on down face
595 //---------------------------
596 el_pt = new TElement<3, 2>;
597
598
599 // LBD
600 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k);
601
602 // RFD
603 el_pt->node_pt(2) = vertex_node_pt(i, j, k);
604
605 // Central down face
607
608 // Central node
609 el_pt->node_pt(3) = central_node_pt(i, j, k);
610
611 Element_pt.push_back(el_pt);
612
613
614 // BACK FACE
615 //===========
616
617 // Left element on back face
618 //---------------------------
619 el_pt = new TElement<3, 2>;
620
621 // LBD
622 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k);
623
624 // Central back face
625 el_pt->node_pt(2) = front_face_node_pt(i, j + 1, k);
626
627 // LBU
628 el_pt->node_pt(1) = vertex_node_pt(i, j + 1, k + 1);
629
630 // central node
631 el_pt->node_pt(3) = central_node_pt(i, j, k);
632
633 Element_pt.push_back(el_pt);
634
635
636 // Down element on back face
637 //---------------------------
638 el_pt = new TElement<3, 2>;
639
640 // LBD
641 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k);
642
643 // RBD
644 el_pt->node_pt(2) = vertex_node_pt(i + 1, j + 1, k);
645
646 // Central back face
647 el_pt->node_pt(1) = front_face_node_pt(i, j + 1, k);
648
649 // Central node
650 el_pt->node_pt(3) = central_node_pt(i, j, k);
651
652 Element_pt.push_back(el_pt);
653
654
655 // Right element on back face
656 //---------------------------
657 el_pt = new TElement<3, 2>;
658
659 // RBD
660 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k);
661
662 // RBU
663 el_pt->node_pt(2) = vertex_node_pt(i + 1, j + 1, k + 1);
664
665 // Central back face
666 el_pt->node_pt(1) = front_face_node_pt(i, j + 1, k);
667
668 // Central node
669 el_pt->node_pt(3) = central_node_pt(i, j, k);
670
671 Element_pt.push_back(el_pt);
672
673
674 // Up element on back face
675 //---------------------------
676 el_pt = new TElement<3, 2>;
677
678 // RBU
679 el_pt->node_pt(0) = vertex_node_pt(i + 1, j + 1, k + 1);
680
681 // LBU
682 el_pt->node_pt(2) = vertex_node_pt(i, j + 1, k + 1);
683
684 // Central back face
685 el_pt->node_pt(1) = front_face_node_pt(i, j + 1, k);
686
687 // Central node
688 el_pt->node_pt(3) = central_node_pt(i, j, k);
689
690 Element_pt.push_back(el_pt);
691
692
693 // LEFT FACE
694 //===========
695
696 // Front element on left face
697 //---------------------------
698 el_pt = new TElement<3, 2>;
699
700 // LFD
701 el_pt->node_pt(0) = vertex_node_pt(i, j, k);
702
703 // Central left face
705
706 // LFU
707 el_pt->node_pt(1) = vertex_node_pt(i, j, k + 1);
708
709 // Central node
710 el_pt->node_pt(3) = central_node_pt(i, j, k);
711
712 Element_pt.push_back(el_pt);
713
714
715 // Down element on left face
716 //---------------------------
717 el_pt = new TElement<3, 2>;
718
719 // LFD
720 el_pt->node_pt(0) = vertex_node_pt(i, j, k);
721
722 // LBD
723 el_pt->node_pt(2) = vertex_node_pt(i, j + 1, k);
724
725 // Central left face
727
728 // central node
729 el_pt->node_pt(3) = central_node_pt(i, j, k);
730
731 Element_pt.push_back(el_pt);
732
733
734 // Back element on left face
735 //---------------------------
736 el_pt = new TElement<3, 2>;
737
738 // LBD
739 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k);
740
741 // LBU
742 el_pt->node_pt(2) = vertex_node_pt(i, j + 1, k + 1);
743
744 // Central left face
746
747 // central node
748 el_pt->node_pt(3) = central_node_pt(i, j, k);
749
750 Element_pt.push_back(el_pt);
751
752
753 // Up element on left face
754 //---------------------------
755 el_pt = new TElement<3, 2>;
756
757 // LBU
758 el_pt->node_pt(0) = vertex_node_pt(i, j + 1, k + 1);
759
760 // LFU
761 el_pt->node_pt(2) = vertex_node_pt(i, j, k + 1);
762
763 // Central left face
765
766 // Central node
767 el_pt->node_pt(3) = central_node_pt(i, j, k);
768
769 Element_pt.push_back(el_pt);
770 }
771 }
772 }
773
774 if (Node_pt.size() != nnod)
775 {
776 std::ostringstream error_stream;
777 error_stream << "Some internal error in the constructor\n"
778 << "Actual number of nodes : " << Node_pt.size()
779 << "\ndoesn't match the prediction : " << nnod
780 << std::endl;
781
782 throw OomphLibError(
784 }
785
786
787 if (Element_pt.size() != nelem)
788 {
789 std::ostringstream error_stream;
790 error_stream << "Some internal error in the constructor\n"
791 << "Actual number of elements : " << Element_pt.size()
792 << "\ndoesn't match the prediction : " << nelem
793 << std::endl;
794
795 throw OomphLibError(
797 }
798
799 // // Deliberately switch nodes to invert elements -- test for self_test()
800 // nelem=nelement();
801 // for (unsigned e=0;e<nelem;e++)
802 // {
803 // Node* node_pt=finite_element_pt(e)->node_pt(0);
804 // finite_element_pt(e)->node_pt(0)=
805 // finite_element_pt(e)->node_pt(1);
806 // finite_element_pt(e)->node_pt(1)=node_pt;
807 // }
808 }
809
810
811} // namespace oomph
cstr elem_len * i
Definition cfortran.h:603
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition mesh.h:183
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition mesh.h:509
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition mesh.h:440
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition mesh.h:186
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
SimpleCubicScaffoldTetMesh(const unsigned &n_x, const unsigned &n_y, const unsigned &n_z, const double &l_x, const double &l_y, const double &l_z, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements and dimensions of cube.
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).