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-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
27#include <algorithm>
28
29#include "map_matrix.h"
30#include "tet_mesh.h"
31#include "Telements.h"
32
33
34namespace oomph
35{
36 //=======================================================================
37 /// Constructor for a FacetedSurface created from a list of nodes
38 /// and connectivity information. This is used in remeshing
39 //=======================================================================
41 Vector<Node*> const& vertex_node_pt,
42 Vector<Vector<unsigned>> const& facet_connectivity,
43 Vector<unsigned> const& facet_boundary_id)
45 {
46 // Create the vertices
47 unsigned n_vertex = vertex_node_pt.size();
48 Vertex_pt.resize(n_vertex);
49 for (unsigned v = 0; v < n_vertex; ++v)
50 {
51 Vertex_pt[v] = new TetMeshVertex(vertex_node_pt[v]);
52 }
53
54 // Create the facets
55 unsigned n_facet = facet_connectivity.size();
56 Facet_pt.resize(n_facet);
57 for (unsigned f = 0; f < n_facet; ++f)
58 {
59 unsigned n_vertex_on_facet = facet_connectivity[f].size();
60 Facet_pt[f] = new TetMeshFacet(n_vertex_on_facet);
61 for (unsigned i = 0; i < n_vertex_on_facet; ++i)
62 {
63 Facet_pt[f]->set_vertex_pt(i, Vertex_pt[facet_connectivity[f][i]]);
64 }
65 // Add in the boundary id
66 Facet_pt[f]->set_one_based_boundary_id(facet_boundary_id[f]);
67 }
68 }
69
70 //=================================================================
71 /// Destructor. Delete allocated memory
72 //================================================================
74 {
75 // Delete the facets and the vertices
76 unsigned n_facet = this->nfacet();
77 for (unsigned f = 0; f < n_facet; f++)
78 {
79 delete Facet_pt[f];
80 }
81 unsigned n_vertex = this->nvertex();
82 for (unsigned v = 0; v < n_vertex; v++)
83 {
84 delete Vertex_pt[v];
85 }
86 }
87
88
89 //================================================================
90 /// Global static data that specifies the permitted
91 /// error in the setup of the boundary coordinates
92 //================================================================
94
95
96 //================================================================
97 /// Setup lookup schemes which establish which elements are located
98 /// next to which boundaries (Doc to outfile if it's open).
99 //================================================================
101 {
102 // Should we document the output here
103 bool doc = false;
104
105 if (outfile) doc = true;
106
107 // Number of boundaries
108 unsigned nbound = nboundary();
109
110 // Wipe/allocate storage for arrays
111 Boundary_element_pt.clear();
113 Boundary_element_pt.resize(nbound);
114 Face_index_at_boundary.resize(nbound);
116
117 // Temporary vector of vectors of pointers to elements on the boundaries:
118 // This is a vector to ensure UNIQUE ordering in all processors
119 Vector<Vector<FiniteElement*>> vector_of_boundary_element_pt;
120 vector_of_boundary_element_pt.resize(nbound);
121
122 // Matrix map for working out the fixed face for elements on boundary
124
125 // Loop over elements
126 //-------------------
127 unsigned nel = nelement();
128
129
130 // Get pointer to vector of boundaries that the
131 // node lives on
132 Vector<std::set<unsigned>*> boundaries_pt(4, 0);
133
134 for (unsigned e = 0; e < nel; e++)
135 {
136 // Get pointer to element
138
139 if (doc) outfile << "Element: " << e << " " << fe_pt << std::endl;
140
141 // Only include 3D elements! Some meshes contain interface elements too.
142 if (fe_pt->dim() == 3)
143 {
144 // Loop over the element's nodes and find out which boundaries they're
145 // on
146 // ----------------------------------------------------------------------
147 // We need only loop over the corner nodes
148 for (unsigned i = 0; i < 4; i++)
149 {
150 fe_pt->node_pt(i)->get_boundaries_pt(boundaries_pt[i]);
151 }
152
153 // Find the common boundaries of each face
155
156 // NOTE: Face indices defined in Telements.h
157
158 // Face 3 connnects points 0, 1 and 2
159 if (boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[2])
160 {
161 std::set<unsigned> aux;
162
163 std::set_intersection(
164 boundaries_pt[0]->begin(),
165 boundaries_pt[0]->end(),
166 boundaries_pt[1]->begin(),
167 boundaries_pt[1]->end(),
168 std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
169
170 std::set_intersection(
171 aux.begin(),
172 aux.end(),
173 boundaries_pt[2]->begin(),
174 boundaries_pt[2]->end(),
175 std::insert_iterator<std::set<unsigned>>(face[3], face[3].begin()));
176 }
177
178 // Face 2 connects points 0, 1 and 3
179 if (boundaries_pt[0] && boundaries_pt[1] && boundaries_pt[3])
180 {
181 std::set<unsigned> aux;
182
183 std::set_intersection(
184 boundaries_pt[0]->begin(),
185 boundaries_pt[0]->end(),
186 boundaries_pt[1]->begin(),
187 boundaries_pt[1]->end(),
188 std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
189
190 std::set_intersection(
191 aux.begin(),
192 aux.end(),
193 boundaries_pt[3]->begin(),
194 boundaries_pt[3]->end(),
195 std::insert_iterator<std::set<unsigned>>(face[2], face[2].begin()));
196 }
197
198 // Face 1 connects points 0, 2 and 3
199 if (boundaries_pt[0] && boundaries_pt[2] && boundaries_pt[3])
200 {
201 std::set<unsigned> aux;
202
203 std::set_intersection(
204 boundaries_pt[0]->begin(),
205 boundaries_pt[0]->end(),
206 boundaries_pt[2]->begin(),
207 boundaries_pt[2]->end(),
208 std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
209
210 std::set_intersection(
211 aux.begin(),
212 aux.end(),
213 boundaries_pt[3]->begin(),
214 boundaries_pt[3]->end(),
215 std::insert_iterator<std::set<unsigned>>(face[1], face[1].begin()));
216 }
217
218 // Face 0 connects points 1, 2 and 3
219 if (boundaries_pt[1] && boundaries_pt[2] && boundaries_pt[3])
220 {
221 std::set<unsigned> aux;
222
223 std::set_intersection(
224 boundaries_pt[1]->begin(),
225 boundaries_pt[1]->end(),
226 boundaries_pt[2]->begin(),
227 boundaries_pt[2]->end(),
228 std::insert_iterator<std::set<unsigned>>(aux, aux.begin()));
229
230 std::set_intersection(
231 aux.begin(),
232 aux.end(),
233 boundaries_pt[3]->begin(),
234 boundaries_pt[3]->end(),
235 std::insert_iterator<std::set<unsigned>>(face[0], face[0].begin()));
236 }
237
238
239 // We now know whether any faces lay on the boundaries
240 for (unsigned i = 0; i < 4; i++)
241 {
242 // How many boundaries are there
243 unsigned count = 0;
244
245 // The number of the boundary
246 int boundary = -1;
247
248 // Loop over all the members of the set and add to the count
249 // and set the boundary
250 for (std::set<unsigned>::iterator it = face[i].begin();
251 it != face[i].end();
252 ++it)
253 {
254 ++count;
255 boundary = *it;
256 }
257
258 // If we're on more than one boundary, this is weird, so die
259 if (count > 1)
260 {
261 std::ostringstream error_stream;
262 fe_pt->output(error_stream);
263 error_stream << "Face " << i << " is on " << count
264 << " boundaries.\n";
265 error_stream << "This is rather strange.\n";
266 error_stream << "Your mesh may be too coarse or your tet mesh\n";
267 error_stream << "may be screwed up. I'm skipping the automated\n";
268 error_stream << "setup of the elements next to the boundaries\n";
269 error_stream << "lookup schemes.\n";
270 OomphLibWarning(error_stream.str(),
271 OOMPH_CURRENT_FUNCTION,
272 OOMPH_EXCEPTION_LOCATION);
273 }
274
275 // If we have a boundary then add this to the appropriate set
276 if (boundary >= 0)
277 {
278 // Does the pointer already exits in the vector
279 Vector<FiniteElement*>::iterator b_el_it = std::find(
280 vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
281 .begin(),
282 vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
283 .end(),
284 fe_pt);
285
286 // Only insert if we have not found it (i.e. got to the end)
287 if (b_el_it ==
288 vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
289 .end())
290 {
291 vector_of_boundary_element_pt[static_cast<unsigned>(boundary)]
292 .push_back(fe_pt);
293 }
294
295 // Also set the fixed face
296 face_identifier(static_cast<unsigned>(boundary), fe_pt) = i;
297 }
298 }
299
300 // Now we set the pointers to the boundary sets to zero
301 for (unsigned i = 0; i < 4; i++)
302 {
303 boundaries_pt[i] = 0;
304 }
305 }
306 }
307
308 // Now copy everything across into permanent arrays
309 //-------------------------------------------------
310
311 // Loop over boundaries
312 //---------------------
313 for (unsigned i = 0; i < nbound; i++)
314 {
315 // Number of elements on this boundary (currently stored in a set)
316 unsigned nel = vector_of_boundary_element_pt[i].size();
317
318 // Allocate storage for the coordinate identifiers
319 Face_index_at_boundary[i].resize(nel);
320
321 unsigned e_count = 0;
323 for (IT it = vector_of_boundary_element_pt[i].begin();
324 it != vector_of_boundary_element_pt[i].end();
325 it++)
326 {
327 // Recover pointer to element
328 FiniteElement* fe_pt = *it;
329
330 // Add to permanent storage
331 Boundary_element_pt[i].push_back(fe_pt);
332
333 Face_index_at_boundary[i][e_count] = face_identifier(i, fe_pt);
334
335 // Increment counter
336 e_count++;
337 }
338 }
339
340
341 // Doc?
342 //-----
343 if (doc)
344 {
345 // Loop over boundaries
346 for (unsigned i = 0; i < nbound; i++)
347 {
348 unsigned nel = Boundary_element_pt[i].size();
349 outfile << "Boundary: " << i << " is adjacent to " << nel << " elements"
350 << std::endl;
351
352 // Loop over elements on given boundary
353 for (unsigned e = 0; e < nel; e++)
354 {
356 outfile << "Boundary element:" << fe_pt
357 << " Face index of boundary is "
358 << Face_index_at_boundary[i][e] << std::endl;
359 }
360 }
361 }
362
363
364 // Lookup scheme has now been setup yet
366 }
367
368
369 //======================================================================
370 /// Assess mesh quality: Ratio of max. edge length to min. height,
371 /// so if it's very large it's BAAAAAD.
372 //======================================================================
373 void TetMeshBase::assess_mesh_quality(std::ofstream& some_file)
374 {
375 Vector<Vector<double>> edge(6);
376 for (unsigned e = 0; e < 6; e++)
377 {
378 edge[e].resize(3);
379 }
380 unsigned nel = this->nelement();
381 for (unsigned e = 0; e < nel; e++)
382 {
383 FiniteElement* fe_pt = this->finite_element_pt(e);
384 for (unsigned i = 0; i < 3; i++)
385 {
386 edge[0][i] = fe_pt->node_pt(2)->x(i) - fe_pt->node_pt(1)->x(i);
387 edge[1][i] = fe_pt->node_pt(0)->x(i) - fe_pt->node_pt(2)->x(i);
388 edge[2][i] = fe_pt->node_pt(1)->x(i) - fe_pt->node_pt(0)->x(i);
389 edge[3][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(0)->x(i);
390 edge[4][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(1)->x(i);
391 edge[5][i] = fe_pt->node_pt(3)->x(i) - fe_pt->node_pt(2)->x(i);
392 }
393
394 double max_length = 0.0;
395 for (unsigned j = 0; j < 6; j++)
396 {
397 double length = 0.0;
398 for (unsigned i = 0; i < 3; i++)
399 {
400 length += edge[j][i] * edge[j][i];
401 }
402 length = sqrt(length);
403 if (length > max_length) max_length = length;
404 }
405
406
407 double min_height = DBL_MAX;
408 for (unsigned j = 0; j < 4; j++)
409 {
410 Vector<double> normal(3);
411 unsigned e0 = 0;
412 unsigned e1 = 0;
413 unsigned e2 = 0;
414 switch (j)
415 {
416 case 0:
417 e0 = 4;
418 e1 = 5;
419 e2 = 1;
420 break;
421
422 case 1:
423 e0 = 1;
424 e1 = 3;
425 e2 = 2;
426 break;
427
428 case 2:
429 e0 = 3;
430 e1 = 4;
431 e2 = 1;
432 break;
433
434 case 3:
435 e0 = 1;
436 e1 = 2;
437 e2 = 3;
438 break;
439
440 default:
441
442 oomph_info << "never get here\n";
443 abort();
444 }
445
446 normal[0] = edge[e0][1] * edge[e1][2] - edge[e0][2] * edge[e1][1];
447 normal[1] = edge[e0][2] * edge[e1][0] - edge[e0][0] * edge[e1][2];
448 normal[2] = edge[e0][0] * edge[e1][1] - edge[e0][1] * edge[e1][0];
449 double norm =
450 normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
451 double inv_norm = 1.0 / sqrt(norm);
452 normal[0] *= inv_norm;
453 normal[1] *= inv_norm;
454 normal[2] *= inv_norm;
455
456 double height = fabs(edge[e2][0] * normal[0] + edge[e2][1] * normal[1] +
457 edge[e2][2] * normal[2]);
458
459 if (height < min_height) min_height = height;
460 }
461
462 double aspect_ratio = max_length / min_height;
463
464 some_file << "ZONE N=4, E=1, F=FEPOINT, ET=TETRAHEDRON\n";
465 for (unsigned j = 0; j < 4; j++)
466 {
467 for (unsigned i = 0; i < 3; i++)
468 {
469 some_file << fe_pt->node_pt(j)->x(i) << " ";
470 }
471 some_file << aspect_ratio << std::endl;
472 }
473 some_file << "1 2 3 4" << std::endl;
474 }
475 }
476
477
478 //======================================================================
479 /// Move the nodes on boundaries with associated Geometric Objects (if any)
480 /// so that they exactly coincide with the geometric object. This requires
481 /// that the boundary coordinates are set up consistently
482 //======================================================================
484 {
485 // Backup in case elements get inverted
486 std::map<Node*, Vector<double>> old_nodal_posn;
487 std::map<Node*, Vector<double>> new_nodal_posn;
488 unsigned nnod = nnode();
489 for (unsigned j = 0; j < nnod; j++)
490 {
491 Node* nod_pt = node_pt(j);
492 Vector<double> x(3);
493 nod_pt->position(x);
494 old_nodal_posn[nod_pt] = x;
495 }
496
497 // Loop over all boundaries
498 unsigned n_bound = this->nboundary();
499 for (unsigned b = 0; b < n_bound; b++)
500 {
501 bool do_it = true;
502
503 // Accumulate reason for not snapping
504 std::stringstream reason;
505 reason << "Can't snap nodes on boundary " << b
506 << " onto geom object because: \n";
507
508 TetMeshFacetedSurface* faceted_surface_pt = 0;
509 std::map<unsigned, TetMeshFacetedSurface*>::iterator it =
511 if (it != Tet_mesh_faceted_surface_pt.end())
512 {
513 faceted_surface_pt = (*it).second;
514 }
515
516 // Facet associated with this boundary?
517 if (faceted_surface_pt == 0)
518 {
519 reason << "-- no facets asssociated with boundary\n";
520 do_it = false;
521 }
522
523 // Get geom object associated with facet
524 GeomObject* geom_obj_pt = 0;
525 if (do_it)
526 {
527 geom_obj_pt = faceted_surface_pt->geom_object_with_boundaries_pt();
528 if (geom_obj_pt == 0)
529 {
530 reason << "-- no geom object associated with boundary\n";
531 do_it = false;
532 }
533 }
534
535 // Triangular facet?
537 {
538 reason << "-- facet has to be triangular and vertex coordinates have\n"
539 << " to have been set up\n";
540 do_it = false;
541 }
542
543 // We need boundary coordinates!
545 {
546 reason << "-- no boundary coordinates were set up\n";
547 do_it = false;
548 }
549
550
551 // Which facet is associated with this boundary?
552 unsigned facet_id_of_boundary = 0;
553 TetMeshFacet* f_pt = 0;
554 if (do_it)
555 {
556 unsigned nf = faceted_surface_pt->nfacet();
557 for (unsigned f = 0; f < nf; f++)
558 {
559 if ((faceted_surface_pt->one_based_facet_boundary_id(f) - 1) == b)
560 {
561 facet_id_of_boundary = f;
562 break;
563 }
564 }
565 f_pt = faceted_surface_pt->facet_pt(facet_id_of_boundary);
566
567
568 // Three vertices?
569 unsigned nv = f_pt->nvertex();
570 if (nv != 3)
571 {
572 reason << "-- number of facet vertices is " << nv
573 << " rather than 3\n";
574 do_it = false;
575 }
576
577 // Have we set up zeta coordinates in geometric object?
578 if ((f_pt->vertex_pt(0)->zeta_in_geom_object().size() != 2) ||
579 (f_pt->vertex_pt(1)->zeta_in_geom_object().size() != 2) ||
580 (f_pt->vertex_pt(2)->zeta_in_geom_object().size() != 2))
581 {
582 reason << "-- no boundary coordinates were set up\n";
583 do_it = false;
584 }
585 }
586
587
588 // Are we ready to go?
589 if (!do_it)
590 {
591 const bool tell_us_why = false;
592 if (tell_us_why)
593 {
594 oomph_info << reason.str() << std::endl;
595 }
596 }
597 else
598 {
599 // Setup area coordinantes in triangular facet
602
605
608
609 double detT = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3);
610
611
612 // Boundary coordinate (cartesian coordinates inside facet)
613 Vector<double> zeta(2);
614
615 // Loop over all nodes on that boundary
616 const unsigned n_boundary_node = this->nboundary_node(b);
617 for (unsigned n = 0; n < n_boundary_node; ++n)
618 {
619 // Get the boundary node and coordinates
620 Node* const nod_pt = this->boundary_node_pt(b, n);
621 nod_pt->get_coordinates_on_boundary(b, zeta);
622
623 // Now we have zeta, the cartesian boundary coordinates
624 // in the (assumed to be triangular) boundary facet; let's
625 // work out the area coordinates
626 // Notation as in
627 // https://en.wikipedia.org/wiki/Barycentric_coordinate_system
628 double s0 =
629 ((y2 - y3) * (zeta[0] - x3) + (x3 - x2) * (zeta[1] - y3)) / detT;
630 double s1 =
631 ((y3 - y1) * (zeta[0] - x3) + (x1 - x3) * (zeta[1] - y3)) / detT;
632 double s2 = 1.0 - s0 - s1;
633
634 Vector<double> zeta_in_geom_obj(2, 0.0);
635 Vector<double> position_from_geom_obj(3, 0.0);
636
637 // Vertex zeta coordinates
638 Vector<double> zeta_0(2);
639 zeta_0 = f_pt->vertex_pt(0)->zeta_in_geom_object();
640
641 Vector<double> zeta_1(2);
642 zeta_1 = f_pt->vertex_pt(1)->zeta_in_geom_object();
643
644 Vector<double> zeta_2(2);
645 zeta_2 = f_pt->vertex_pt(2)->zeta_in_geom_object();
646
647
648#ifdef PARANOID
649
650 // Compute zeta values of the vertices from parametrisation of
651 // boundaries
652 double tol = 1.0e-12;
653 Vector<double> zeta_from_boundary(2);
654 Vector<double> zeta_vertex(2);
655 for (unsigned v = 0; v < 3; v++)
656 {
657 zeta_vertex = f_pt->vertex_pt(v)->zeta_in_geom_object();
658 for (unsigned alt = 0; alt < 2; alt++)
659 {
660 switch (v)
661 {
662 case 0:
663
664 if (alt == 0)
665 {
666 faceted_surface_pt->boundary_zeta01(
667 facet_id_of_boundary, 0.0, zeta_from_boundary);
668 }
669 else
670 {
671 faceted_surface_pt->boundary_zeta20(
672 facet_id_of_boundary, 1.0, zeta_from_boundary);
673 }
674 break;
675
676 case 1:
677
678 if (alt == 0)
679 {
680 faceted_surface_pt->boundary_zeta01(
681 facet_id_of_boundary, 1.0, zeta_from_boundary);
682 }
683 else
684 {
685 faceted_surface_pt->boundary_zeta12(
686 facet_id_of_boundary, 0.0, zeta_from_boundary);
687 }
688 break;
689
690 case 2:
691
692 if (alt == 0)
693 {
694 faceted_surface_pt->boundary_zeta12(
695 facet_id_of_boundary, 1.0, zeta_from_boundary);
696 }
697 else
698 {
699 faceted_surface_pt->boundary_zeta20(
700 facet_id_of_boundary, 0.0, zeta_from_boundary);
701 }
702 break;
703 }
704
705 double error =
706 sqrt(pow((zeta_vertex[0] - zeta_from_boundary[0]), 2) +
707 pow((zeta_vertex[1] - zeta_from_boundary[1]), 2));
708 if (error > tol)
709 {
710 std::ostringstream error_message;
711 error_message
712 << "Error in parametrisation of boundary coordinates \n"
713 << "for vertex " << v << " [alt=" << alt << "] in facet "
714 << facet_id_of_boundary << " : \n"
715 << "zeta_vertex = [ " << zeta_vertex[0] << " "
716 << zeta_vertex[1] << " ] \n"
717 << "zeta_from_boundary = [ " << zeta_from_boundary[0]
718 << " " << zeta_from_boundary[1] << " ] \n"
719 << std::endl;
720 throw OomphLibError(error_message.str(),
721 OOMPH_CURRENT_FUNCTION,
722 OOMPH_EXCEPTION_LOCATION);
723 }
724 }
725 }
726
727#endif
728
729 // Compute zeta values of the interpolation parameters
730 Vector<double> zeta_a(3, 0.0);
731 Vector<double> zeta_b(3, 0.0);
732 Vector<double> zeta_c(3, 0.0);
733 Vector<double> zeta_d(3, 0.0);
734 Vector<double> zeta_e(3, 0.0);
735 Vector<double> zeta_f(3, 0.0);
736 faceted_surface_pt->boundary_zeta01(facet_id_of_boundary, s1, zeta_a);
737 faceted_surface_pt->boundary_zeta01(
738 facet_id_of_boundary, 1.0 - s0, zeta_d);
739
740 faceted_surface_pt->boundary_zeta12(facet_id_of_boundary, s2, zeta_c);
741 faceted_surface_pt->boundary_zeta12(
742 facet_id_of_boundary, 1.0 - s1, zeta_f);
743
744 faceted_surface_pt->boundary_zeta20(
745 facet_id_of_boundary, 1.0 - s2, zeta_b);
746 faceted_surface_pt->boundary_zeta20(facet_id_of_boundary, s0, zeta_e);
747
748 // Transfinite mapping
749 zeta_in_geom_obj[0] = s0 * (zeta_a[0] + zeta_b[0] - zeta_0[0]) +
750 s1 * (zeta_c[0] + zeta_d[0] - zeta_1[0]) +
751 s2 * (zeta_e[0] + zeta_f[0] - zeta_2[0]);
752 zeta_in_geom_obj[1] = s0 * (zeta_a[1] + zeta_b[1] - zeta_0[1]) +
753 s1 * (zeta_c[1] + zeta_d[1] - zeta_1[1]) +
754 s2 * (zeta_e[1] + zeta_f[1] - zeta_2[1]);
755
756 unsigned n_tvalues =
757 1 + nod_pt->position_time_stepper_pt()->nprev_values();
758 for (unsigned t = 0; t < n_tvalues; ++t)
759 {
760 // Get the position according to the underlying geometric object
761 geom_obj_pt->position(t, zeta_in_geom_obj, position_from_geom_obj);
762
763 // Move the node
764 for (unsigned i = 0; i < 3; i++)
765 {
766 nod_pt->x(t, i) = position_from_geom_obj[i];
767 }
768 }
769 }
770 }
771 }
772
773 // Check if any element is inverted
774 bool some_element_is_inverted = false;
775 unsigned count = 0;
776 unsigned nel = nelement();
777 for (unsigned e = 0; e < nel; e++)
778 {
780 bool passed = true;
781 el_pt->check_J_eulerian_at_knots(passed);
782 if (!passed)
783 {
784 some_element_is_inverted = true;
785 char filename[100];
786 std::ofstream some_file;
787 sprintf(filename, "overly_distorted_element%i.dat", count);
788 some_file.open(filename);
789 unsigned nnod_1d = el_pt->nnode_1d();
790 el_pt->output(some_file, nnod_1d);
791 some_file.close();
792
793 // Reset to old nodal position
794 unsigned nnod = el_pt->nnode();
795 for (unsigned j = 0; j < nnod; j++)
796 {
797 Node* nod_pt = el_pt->node_pt(j);
798 Vector<double> x_current(3);
799 nod_pt->position(x_current);
800 new_nodal_posn[nod_pt] = x_current;
801 Vector<double> old_x(old_nodal_posn[nod_pt]);
802 for (unsigned i = 0; i < 3; i++)
803 {
804 nod_pt->x(i) = old_x[i];
805 }
806 }
807
808 // Plot
809 sprintf(filename, "orig_overly_distorted_element%i.dat", count);
810 some_file.open(filename);
811 el_pt->output(some_file, nnod_1d);
812 some_file.close();
813
814 // Reset
815 for (unsigned j = 0; j < nnod; j++)
816 {
817 Node* nod_pt = el_pt->node_pt(j);
818 for (unsigned i = 0; i < 3; i++)
819 {
820 nod_pt->x(i) = new_nodal_posn[nod_pt][i];
821 }
822 }
823
824 // Bump
825 count++;
826 }
827 }
828 if (some_element_is_inverted)
829 {
830 std::ostringstream error_message;
831 error_message
832 << "A number of elements, namely: " << count
833 << " are inverted after snapping. Their shapes are in "
834 << " overly_distorted_element*.dat and "
835 "orig_overly_distorted_element*.dat"
836 << "Next person to get this error: Please implement a straightforward\n"
837 << "variant of one of the functors in src/mesh_smoothing to switch\n"
838 << "to harmonic mapping\n"
839 << std::endl;
840 throw OomphLibError(
841 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
842 }
843 else
844 {
845 oomph_info << "No elements are inverted after snapping. Yay!"
846 << std::endl;
847 }
848 }
849
850
851 /// //////////////////////////////////////////////////////////////////
852 /// //////////////////////////////////////////////////////////////////
853 /// //////////////////////////////////////////////////////////////////
854
855
856} // namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A general Finite Element class.
Definition: elements.h:1313
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void check_J_eulerian_at_knots(bool &passed) const
Check that Jacobian of mapping between local and Eulerian coordinates at all integration points is po...
Definition: elements.cc:4237
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
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition: mesh.h:91
bool Lookup_for_elements_next_boundary_is_setup
Flag to indicate that the lookup schemes for elements that are adjacent to the boundaries has been se...
Definition: mesh.h:87
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:190
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:95
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
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
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2499
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b....
Definition: nodes.cc:2379
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....
void snap_nodes_onto_geometric_objects()
Move the nodes on boundaries with associated GeomObjects so that they exactly coincide with the geome...
Definition: tet_mesh.cc:483
std::map< unsigned, TetMeshFacetedSurface * > Tet_mesh_faceted_surface_pt
Reverse lookup scheme: Pointer to faceted surface (if any!) associated with boundary b.
Definition: tet_mesh.h:1045
static double Tolerance_for_boundary_finding
Global static data that specifies the permitted error in the setup of the boundary coordinates.
Definition: tet_mesh.h:677
void assess_mesh_quality(std::ofstream &some_file)
Assess mesh quality: Ratio of max. edge length to min. height, so if it's very large it's BAAAAAD.
Definition: tet_mesh.cc:373
void setup_boundary_element_info()
Setup lookup schemes which establish which elements are located next to mesh's boundaries (wrapper to...
Definition: tet_mesh.h:1007
std::map< unsigned, Vector< Vector< double > > > Triangular_facet_vertex_boundary_coordinate
Boundary coordinates of vertices in triangular facets associated with given boundary....
Definition: tet_mesh.h:1054
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:168
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:181
TetMeshVertex * vertex_pt(const unsigned &j) const
Pointer to j-th vertex (const)
Definition: tet_mesh.h:187
virtual ~TetMeshFacetedClosedSurfaceForRemesh()
Destructor. Delete allocated memory.
Definition: tet_mesh.cc:73
TetMeshFacetedClosedSurfaceForRemesh(Vector< Node * > const &vertex_node_pt, Vector< Vector< unsigned > > const &facet_connectivity, Vector< unsigned > const &facet_boundary_id)
Constructor for a FacetedSurface created from a list of nodes and connectivity information....
Definition: tet_mesh.cc:40
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:538
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:306
TetMeshFacet * facet_pt(const unsigned &j) const
Pointer to j-th facet.
Definition: tet_mesh.h:373
virtual void boundary_zeta01(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:417
Vector< TetMeshVertex * > Vertex_pt
Vector pointers to vertices.
Definition: tet_mesh.h:483
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
virtual void boundary_zeta12(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:436
Vector< TetMeshFacet * > Facet_pt
Vector of pointers to facets.
Definition: tet_mesh.h:486
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
DiskLikeGeomObjectWithBoundaries * geom_object_with_boundaries_pt()
Access to GeomObject with boundaries associated with this surface (Null if there isn't one!...
Definition: tet_mesh.h:386
virtual void boundary_zeta20(const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)
Virtual function that specifies the variation of the zeta coordinates in the GeomObject along the bou...
Definition: tet_mesh.h:455
unsigned nvertex() const
Number of vertices.
Definition: tet_mesh.h:319
Vertex for Tet mesh generation. Can lie on multiple boundaries (identified via one-based enumeration!...
Definition: tet_mesh.h:50
Vector< double > zeta_in_geom_object() const
Get intrinisic coordinates in GeomObject (returns zero sized vector if no such coordinates have been ...
Definition: tet_mesh.h:118
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...