mesh.h
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// Header file for general mesh classes
27
28// Include guards to prevent multiple inclusion of the header
29#ifndef OOMPH_GENERIC_MESH_HEADER
30#define OOMPH_GENERIC_MESH_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37#ifdef OOMPH_HAS_MPI
38#include "mpi.h"
39#endif
40
41#include <float.h>
42#include <list>
43#include <typeinfo>
44#include <string>
45
46// oomph-lib headers
47#include "Vector.h"
48#include "nodes.h"
49#include "elements.h"
50#include "timesteppers.h"
52#include "matrices.h"
53#include "refineable_elements.h"
54
55namespace oomph
56{
57 //=================================================================
58 /// A general mesh class.
59 ///
60 /// The main components of a Mesh are:
61 /// - pointers to its Nodes
62 /// - pointers to its Elements
63 /// - pointers to its boundary Nodes
64 ///
65 //=================================================================
66 class Mesh
67 {
68 public:
69 /// Problem is a friend
70 friend class Problem;
71
72
73 /// Default Steady Timestepper, to be used in default arguments
74 /// to Mesh constructors
76
77
78 protected:
79 /// Vector of Vector of pointers to nodes on the boundaries:
80 /// Boundary_node_pt(b,n). Note that this is private to force
81 /// the use of the add_boundary_node() function, which ensures
82 /// that the reverse look-up schemes for the nodes are set up.
84
85 /// Flag to indicate that the lookup schemes for elements that are adjacent
86 /// to the boundaries has been set up.
88
89 /// Vector of Vector of pointers to elements on the boundaries:
90 /// Boundary_element_pt(b,e)
92
93 /// For the e-th finite element on boundary b, this is the index of
94 /// the face that lies along that boundary
96
97#ifdef OOMPH_HAS_MPI
98
99 /// Map of vectors holding the pointers to the root halo elements
100 std::map<unsigned, Vector<GeneralisedElement*>> Root_halo_element_pt;
101
102 /// Map of vectors holding the pointers to the root haloed elements
103 std::map<unsigned, Vector<GeneralisedElement*>> Root_haloed_element_pt;
104
105 /// Map of vectors holding the pointers to the halo nodes
106 std::map<unsigned, Vector<Node*>> Halo_node_pt;
107
108 /// Map of vectors holding the pointers to the haloed nodes
109 std::map<unsigned, Vector<Node*>> Haloed_node_pt;
110
111 /// Map of vectors holding the pointers to the shared nodes.
112 /// These are all the nodes that are on two "neighbouring" processes
113 /// (the halo(ed) lookup scheme depends upon which processor is in charge
114 /// - a node which is on 3 processors, for example, will not feature in
115 /// the halo(ed) lookup scheme between the two lowest-numbered processors)
116 std::map<unsigned, Vector<Node*>> Shared_node_pt;
117
118 /// Pointer to communicator -- set to NULL if mesh is not distributed
120
121 /// External halo(ed) elements are created as and when they are needed
122 /// to act as source elements for the particular process's mesh.
123 /// The storage is wiped and rebuilt every time the mesh is refined.
124
125 /// Map of vectors holding the pointers to the external halo elements
126 std::map<unsigned, Vector<GeneralisedElement*>> External_halo_element_pt;
127
128 /// Map of vectors holding the pointers to the external haloed elements
129 std::map<unsigned, Vector<GeneralisedElement*>> External_haloed_element_pt;
130
131
132 // External halo(ed) nodes are on the external halo(ed) elements
133
134 /// Map of vectors holding the pointers to the external halo nodes
135 std::map<unsigned, Vector<Node*>> External_halo_node_pt;
136
137 /// Map of vectors holding the pointers to the external haloed nodes
138 std::map<unsigned, Vector<Node*>> External_haloed_node_pt;
139
140 /// bool to indicate whether to keep all elements in a mesh as halos or not
142
143 /// Set this to true to suppress resizing of halo nodes (at your own risk!)
145
146 /// Setup shared node scheme
148
149#endif
150
151 /// Assign the global equation numbers in the Data stored at the
152 /// nodes and also internal element Data. Also, build (via push_back) the
153 /// Vector of pointers to the dofs (variables).
154 unsigned long assign_global_eqn_numbers(Vector<double*>& Dof_pt);
155
156 /// Function to describe the dofs of the Mesh. The ostream
157 /// specifies the output stream to which the description
158 /// is written; the string stores the currently
159 /// assembled output that is ultimately written to the
160 /// output stream by Data::describe_dofs(...); it is typically
161 /// built up incrementally as we descend through the
162 /// call hierarchy of this function when called from
163 /// Problem::describe_dofs(...)
164 void describe_dofs(std::ostream& out,
165 const std::string& current_string) const;
166
167 /// Function to describe the local dofs of the elements. The ostream
168 /// specifies the output stream to which the description
169 /// is written; the string stores the currently
170 /// assembled output that is ultimately written to the
171 /// output stream by Data::describe_dofs(...); it is typically
172 /// built up incrementally as we descend through the
173 /// call hierarchy of this function when called from
174 /// Problem::describe_dofs(...)
175 void describe_local_dofs(std::ostream& out,
176 const std::string& current_string) const;
177
178 /// Assign the local equation numbers in all elements
179 /// If the boolean argument is true then also store pointers to dofs
180 void assign_local_eqn_numbers(const bool& store_local_dof_pt);
181
182 /// Vector of pointers to nodes
184
185 /// Vector of pointers to generalised elements
187
188 /// Vector of boolean data that indicates whether the boundary
189 /// coordinates have been set for the boundary
190 std::vector<bool> Boundary_coordinate_exists;
191
192 /// A function that upgrades an ordinary node to a boundary node
193 /// We shouldn't ever really use this, but it does make life that
194 /// bit easier for the lazy mesh writer. The pointer to the node is
195 /// replaced by a pointer to the new boundary node in all element look-up
196 /// schemes and in the mesh's Node_pt vector. The new node is also
197 /// addressed by node_pt on return from the function.
200
202
203
204 public:
205#ifdef OOMPH_HAS_MPI
206
207
208 /// Helper function that resizes halo nodes to the same
209 /// size as their non-halo counterparts if required. (A discrepancy
210 /// can arise if a FaceElement that introduces additional unknowns
211 /// are attached to a bulk element that shares a node with a haloed element.
212 /// In that case the joint node between haloed and non-haloed element
213 /// is resized on that processor but not on the one that holds the
214 /// halo counterpart (because no FaceElement is attached to the halo
215 /// element)
216 void resize_halo_nodes();
217
218#endif
219
220
221 /// Typedef for function pointer to function that computes
222 /// steady exact solution
224 const Vector<double>& x, Vector<double>& soln);
225
226 /// Typedef for function pointer to function that computes unsteady
227 /// exact solution
229 const double& time, const Vector<double>& x, Vector<double>& soln);
230
231 /// Boolean used to control warning about empty mesh level
232 /// timestepper function
234
235 /// Default constructor
237 {
238 // Lookup scheme hasn't been setup yet
240#ifdef OOMPH_HAS_MPI
241 // Set defaults for distributed meshes
242
243 // Mesh hasn't been distributed: Null out pointer to communicator
244 Comm_pt = 0;
245 // Don't keep all objects as halos
247 // Don't output halo elements
248 Output_halo_elements = false;
249 // Don't suppress automatic resizing of halo nodes
251#endif
252 }
253
254 /// Constructor builds combined mesh from the meshes specified.
255 /// Note: This simply merges the meshes' elements and nodes (ignoring
256 /// duplicates; no boundary information etc. is created).
257 Mesh(const Vector<Mesh*>& sub_mesh_pt)
258 {
259#ifdef OOMPH_HAS_MPI
260 // Mesh hasn't been distributed: Null out pointer to communicator
261 Comm_pt = 0;
262#endif
263 // Now merge the meshes
264 merge_meshes(sub_mesh_pt);
265 }
266
267 /// Merge meshes.
268 /// Note: This simply merges the meshes' elements and nodes (ignoring
269 /// duplicates; no boundary information etc. is created).
270 void merge_meshes(const Vector<Mesh*>& sub_mesh_pt);
271
272 /// Interface for function that is used to setup the boundary
273 /// information (Empty virtual function -- implement this for specific
274 /// Mesh classes)
276
277 /// Setup lookup schemes which establish whic elements are located
278 /// next to mesh's boundaries. Doc in outfile (if it's open).
279 /// (Empty virtual function -- implement this for specific
280 /// Mesh classes)
281 virtual void setup_boundary_element_info(std::ostream& outfile) {}
282
283 /// Virtual function to perform the reset boundary elements info rutines
285 Vector<unsigned>& ntmp_boundary_elements,
286 Vector<Vector<unsigned>>& ntmp_boundary_elements_in_region,
287 Vector<FiniteElement*>& deleted_elements)
288 {
289 std::ostringstream error_stream;
290 error_stream << "Empty default reset boundary element info function"
291 << "called.\n";
292 error_stream << "This should be overloaded in a specific "
293 << "TriangleMeshBase\n";
294 throw OomphLibError(error_stream.str(),
295 "Mesh::reset_boundary_element_info()",
296 OOMPH_EXCEPTION_LOCATION);
297 }
298
299 /// Output boundary coordinates on boundary b -- template argument
300 /// specifies the bulk element type (needed to create FaceElement
301 /// of appropriate type on mesh boundary).
302 template<class BULK_ELEMENT>
303 void doc_boundary_coordinates(const unsigned& b, std::ofstream& the_file)
304 {
305 if (nelement() == 0) return;
307 {
308 oomph_info << "No boundary coordinates were set up for boundary " << b
309 << std::endl;
310 return;
311 }
312
313 // Get spatial dimension
314 unsigned dim = finite_element_pt(0)->node_pt(0)->ndim();
315
316 // Loop over all elements on boundaries
317 unsigned nel = this->nboundary_element(b);
318
319 // Loop over the bulk elements adjacent to boundary b
320 for (unsigned e = 0; e < nel; e++)
321 {
322 // Get pointer to the bulk element that is adjacent to boundary b
323 FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
324
325 // Find the index of the face of element e along boundary b
326 int face_index = this->face_index_at_boundary(b, e);
327
328 // Create new face element
330 new DummyFaceElement<BULK_ELEMENT>(bulk_elem_pt, face_index);
331
332 // Specify boundary id in bulk mesh (needed to extract
333 // boundary coordinate)
335
336 // Doc boundary coordinate
337 Vector<double> s(dim - 1);
338 Vector<double> zeta(dim - 1);
339 Vector<double> x(dim);
340 unsigned n_plot = 5;
341 the_file << el_pt->tecplot_zone_string(n_plot);
342
343 // Loop over plot points
344 unsigned num_plot_points = el_pt->nplot_points(n_plot);
345 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
346 {
347 // Get local coordinates of plot point
348 el_pt->get_s_plot(iplot, n_plot, s);
349 el_pt->interpolated_zeta(s, zeta);
350 el_pt->interpolated_x(s, x);
351 for (unsigned i = 0; i < dim; i++)
352 {
353 the_file << x[i] << " ";
354 }
355 for (unsigned i = 0; i < (dim - 1); i++)
356 {
357 the_file << zeta[i] << " ";
358 }
359
360 the_file << std::endl;
361 }
362 el_pt->write_tecplot_zone_footer(the_file, n_plot);
363
364 // Cleanup
365 delete el_pt;
366 }
367 }
368
369
370 /// Scale all nodal coordinates by given factor. Virtual
371 /// so it can be overloaded in SolidMesh class where it also
372 /// re-assigns the Lagrangian coordinates.
373 virtual void scale_mesh(const double& factor)
374 {
375 unsigned nnod = this->nnode();
376 unsigned dim = this->node_pt(0)->ndim();
377 for (unsigned j = 0; j < nnod; j++)
378 {
379 Node* nod_pt = this->node_pt(j);
380 for (unsigned i = 0; i < dim; i++)
381 {
382 nod_pt->x(i) *= factor;
383 }
384 }
385 }
386
387 /// Broken copy constructor
388 Mesh(const Mesh& dummy) = delete;
389
390 /// Broken assignment operator
391 void operator=(const Mesh&) = delete;
392
393 /// Virtual Destructor to clean up all memory
394 virtual ~Mesh();
395
396
397 /// Flush storage for elements and nodes by emptying the
398 /// vectors that store the pointers to them. This is
399 /// useful if a particular mesh is only built to generate
400 /// a small part of a bigger mesh. Once the elements and
401 /// nodes have been created, they are typically copied
402 /// into the new mesh and the auxiliary mesh can be
403 /// deleted. However, if we simply call the destructor
404 /// of the auxiliary mesh, it will also wipe out
405 /// the nodes and elements, because it still "thinks"
406 /// it's in charge of these...
408 {
411 }
412
413 /// Flush storage for elements (only) by emptying the
414 /// vectors that store the pointers to them. This is
415 /// useful if a particular mesh is only built to generate
416 /// a small part of a bigger mesh. Once the elements and
417 /// nodes have been created, they are typically copied
418 /// into the new mesh and the auxiliary mesh can be
419 /// deleted. However, if we simply call the destructor
420 /// of the auxiliary mesh, it will also wipe out
421 /// the nodes and elements, because it still "thinks"
422 /// it's in charge of these...
424 {
425 Element_pt.clear();
426 }
427
428 /// Flush storage for nodes (only) by emptying the
429 /// vectors that store the pointers to them.
431 {
432 Node_pt.clear();
433 }
434
435 /// Return pointer to global node n
436 Node*& node_pt(const unsigned long& n)
437 {
438 return Node_pt[n];
439 }
440
441 /// Return pointer to global node n (const version)
442 Node* node_pt(const unsigned long& n) const
443 {
444 return Node_pt[n];
445 }
446
447 /// Return pointer to element e
448 GeneralisedElement*& element_pt(const unsigned long& e)
449 {
450 return Element_pt[e];
451 }
452
453 /// Return pointer to element e (const version)
454 GeneralisedElement* element_pt(const unsigned long& e) const
455 {
456 return Element_pt[e];
457 }
458
459 /// Return reference to the Vector of elements
461 {
462 return Element_pt;
463 }
464
465 /// Return reference to the Vector of elements
467 {
468 return Element_pt;
469 }
470
471 /// Upcast (downcast?) to FiniteElement
472 /// (needed to access FiniteElement member functions).
473 FiniteElement* finite_element_pt(const unsigned& e) const
474 {
475#ifdef PARANOID
476 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
477 if (el_pt == 0)
478 {
479 // Error
480 throw OomphLibError("Failed cast to FiniteElement* ",
481 OOMPH_CURRENT_FUNCTION,
482 OOMPH_EXCEPTION_LOCATION);
483 // Dummy return to keep intel compiler happy
484 return el_pt;
485 }
486 return el_pt;
487#else
488 return dynamic_cast<FiniteElement*>(Element_pt[e]);
489#endif
490 }
491
492 /// Return pointer to node n on boundary b
493 Node*& boundary_node_pt(const unsigned& b, const unsigned& n)
494 {
495 return Boundary_node_pt[b][n];
496 }
497
498 /// Return pointer to node n on boundary b
499 Node* boundary_node_pt(const unsigned& b, const unsigned& n) const
500 {
501 return Boundary_node_pt[b][n];
502 }
503
504 /// Set the number of boundaries in the mesh
505 void set_nboundary(const unsigned& nbound)
506 {
507 Boundary_node_pt.resize(nbound);
508 Boundary_coordinate_exists.resize(nbound, false);
509 }
510
511 /// Clear all pointers to boundary nodes
513
514 /// Remove all information about nodes stored on the b-th
515 /// boundary of the mesh
516 void remove_boundary_nodes(const unsigned& b);
517
518 /// Remove a node from the boundary b
519 void remove_boundary_node(const unsigned& b, Node* const& node_pt);
520
521 /// Add a (pointer to) a node to the b-th boundary
522 void add_boundary_node(const unsigned& b, Node* const& node_pt);
523
524 /// Replace existing boundary node lookup schemes with new schemes created
525 /// using the boundary data stored in the nodes.
527 {
528 // Clear existing boundary data
529 Boundary_node_pt.clear();
530
531 // Loop over nodes adding them to the appropriate boundary lookup schemes
532 // in the mesh.
533 const unsigned n_node = nnode();
534 for (unsigned nd = 0; nd < n_node; nd++)
535 {
536 Node* nd_pt = node_pt(nd);
537
538 if (nd_pt->is_on_boundary())
539 {
540 // Get set of boundaries that the node is on
541 std::set<unsigned>* boundaries_pt;
542 nd_pt->get_boundaries_pt(boundaries_pt);
543
544 // If needed then add more boundaries to this mesh
545 unsigned max_boundary_n =
546 1 + *std::max_element(boundaries_pt->begin(), boundaries_pt->end());
547 if (max_boundary_n > nboundary())
548 {
549 set_nboundary(max_boundary_n);
550 }
551
552 // Add node pointer to the appropriate Boundary_node_pt vectors
553 std::set<unsigned>::const_iterator it;
554 for (it = boundaries_pt->begin(); it != boundaries_pt->end(); it++)
555 {
556 Boundary_node_pt[*it].push_back(nd_pt);
557 }
558 }
559 }
560 }
561
562 /// Indicate whether the i-th boundary has an intrinsic coordinate.
563 // By default, if the array Boundary_coordinate has not been resized,
564 // return false
565 bool boundary_coordinate_exists(const unsigned& i) const
566 {
567 if (Boundary_coordinate_exists.empty())
568 {
569 return false;
570 }
571 // ALH: This bounds-checking code needs to remain, because
572 // Boundary_coordinate_exists is
573 // an stl vector not our overloaded Vector class
574#ifdef RANGE_CHECKING
575 if (i >= Boundary_coordinate_exists.size())
576 {
577 std::ostringstream error_message;
578 error_message << "Range Error: " << i << " is not in the range (0,"
579 << Boundary_coordinate_exists.size() - 1 << std::endl;
580
581 throw OomphLibError(error_message.str(),
582 OOMPH_CURRENT_FUNCTION,
583 OOMPH_EXCEPTION_LOCATION);
584 }
585#endif
587 }
588
589 /// Return number of elements in the mesh
590 unsigned long nelement() const
591 {
592 return Element_pt.size();
593 }
594
595 /// Return number of nodes in the mesh
596 unsigned long nnode() const
597 {
598 return Node_pt.size();
599 }
600
601 /// Return number of dof types in mesh
602 unsigned ndof_types() const;
603
604 /// Return number of elemental dimension in mesh
605 unsigned elemental_dimension() const;
606
607 /// Return number of nodal dimension in mesh
608 unsigned nodal_dimension() const;
609
610 /// Add a (pointer to a) node to the mesh
612 {
613 Node_pt.push_back(node_pt);
614 }
615
616 /// Add a (pointer to) an element to the mesh
618 {
619 Element_pt.push_back(element_pt);
620 }
621
622 /// Update nodal positions in response to changes in the domain
623 /// shape. Uses the FiniteElement::get_x(...) function for FiniteElements
624 /// and doesn't do anything for other element types.
625 /// If a MacroElement pointer has been set for a FiniteElement,
626 /// the MacroElement representation is used to update the
627 /// nodal positions; if not get_x(...) uses the FE interpolation
628 /// and thus leaves the nodal positions unchanged.
629 /// Virtual, so it can be overloaded by specific meshes,
630 /// such as AlgebraicMeshes or SpineMeshes.
631 /// Generally, this function updates the position of all nodes
632 /// in response to changes in the boundary position.
633 /// However, we ignore all SolidNodes since their
634 /// position is computed as part of the solution -- unless
635 /// the bool flag is set to true. Such calls are typically made
636 /// when the initial mesh is created and/or after a mesh has been
637 /// refined repeatedly before the start of the computation.
638 virtual void node_update(const bool& update_all_solid_nodes = false);
639
640 /// Re-order nodes in the order in which they appear in elements --
641 /// can be overloaded for more efficient re-ordering
642 virtual void reorder_nodes(const bool& use_old_ordering = true);
643
644 /// Get a reordering of the nodes in the order in which they
645 /// appear in elements -- can be overloaded for more efficient
646 /// re-ordering
647 virtual void get_node_reordering(Vector<Node*>& reordering,
648 const bool& use_old_ordering = true) const;
649
650 /// Constuct a Mesh of FACE_ELEMENTs along the b-th boundary
651 /// of the mesh (which contains elements of type BULK_ELEMENT)
652 template<class BULK_ELEMENT, template<class> class FACE_ELEMENT>
653 void build_face_mesh(const unsigned& b, Mesh* const& face_mesh_pt)
654 {
655 // Find the number of nodes on the boundary
656 unsigned nbound_node = nboundary_node(b);
657 // Loop over the boundary nodes and add them to face mesh node pointer
658 for (unsigned n = 0; n < nbound_node; n++)
659 {
660 face_mesh_pt->add_node_pt(boundary_node_pt(b, n));
661 }
662
663 // Find the number of elements next to the boundary
664 unsigned nbound_element = nboundary_element(b);
665 // Loop over the elements adjacent to boundary b
666 for (unsigned e = 0; e < nbound_element; e++)
667 {
668 // Create the FaceElement
669 FACE_ELEMENT<BULK_ELEMENT>* face_element_pt =
670 new FACE_ELEMENT<BULK_ELEMENT>(boundary_element_pt(b, e),
672
673 // Add the face element to the face mesh
674 face_mesh_pt->add_element_pt(face_element_pt);
675 }
676
677#ifdef OOMPH_HAS_MPI
678 // If the bulk mesh has been distributed then the face mesh is too
679 if (this->is_mesh_distributed())
680 {
681 face_mesh_pt->set_communicator_pt(this->communicator_pt());
682 }
683#endif
684 }
685
686 /// Self-test: Check elements and nodes. Return 0 for OK
687 unsigned self_test();
688
689
690 /// Determine max and min area for all FiniteElements in the mesh
691 /// (non-FiniteElements are ignored)
692 void max_and_min_element_size(double& max_size, double& min_size)
693 {
694 max_size = 0.0;
695 min_size = DBL_MAX;
696 unsigned nel = nelement();
697 for (unsigned e = 0; e < nel; e++)
698 {
699 max_size = std::max(max_size, finite_element_pt(e)->size());
700 min_size = std::min(min_size, finite_element_pt(e)->size());
701 }
702 }
703
704
705 /// Determine the sum of all "sizes" of the FiniteElements in the
706 /// mesh (non-FiniteElements are ignored). This gives the length/area/volume
707 /// occupied by the mesh
708 double total_size()
709 {
710 double size = 0.0;
711 unsigned nel = nelement();
712 for (unsigned e = 0; e < nel; e++)
713 {
715 if (fe_pt != 0)
716 {
717 size += fe_pt->size();
718 }
719 }
720 return size;
721 }
722
723
724 /// Check for inverted elements and report outcome
725 /// in boolean variable. This visits all elements at their
726 /// integration points and checks if the Jacobian of the
727 /// mapping between local and global coordinates is positive --
728 /// using the same test that would be carried out (but only in PARANOID
729 /// mode) during the assembly of the elements' Jacobian matrices.
730 /// Inverted elements are output in inverted_element_file (if the
731 /// stream is open).
732 void check_inverted_elements(bool& mesh_has_inverted_elements,
733 std::ofstream& inverted_element_file);
734
735
736 /// Check for inverted elements and report outcome
737 /// in boolean variable. This visits all elements at their
738 /// integration points and checks if the Jacobian of the
739 /// mapping between local and global coordinates is positive --
740 /// using the same test that would be carried out (but only in PARANOID
741 /// mode) during the assembly of the elements' Jacobian matrices.
742 void check_inverted_elements(bool& mesh_has_inverted_elements)
743 {
744 std::ofstream inverted_element_file;
745 check_inverted_elements(mesh_has_inverted_elements,
746 inverted_element_file);
747 }
748
749
750 /// Check for repeated nodes within a given spatial tolerance.
751 /// Return (0/1) for (pass/fail).
752 unsigned check_for_repeated_nodes(const double& epsilon = 1.0e-12)
753 {
754 oomph_info << "\n\nStarting check for repeated nodes...";
755 bool failed = false;
756 unsigned nnod = nnode();
757 for (unsigned j = 0; j < nnod; j++)
758 {
759 Node* nod1_pt = this->node_pt(j);
760 unsigned dim = nod1_pt->ndim();
761 for (unsigned k = j + 1; k < nnod; k++)
762 {
763 Node* nod2_pt = this->node_pt(k);
764 double dist = 0.0;
765 for (unsigned i = 0; i < dim; i++)
766 {
767 dist += pow((nod1_pt->x(i) - nod2_pt->x(i)), 2);
768 }
769 dist = sqrt(dist);
770 if (dist < epsilon)
771 {
772 oomph_info << "\n\nRepeated node!" << std::endl;
773 oomph_info << "Distance between nodes " << j << " and " << k
774 << std::endl;
775 oomph_info << "is " << dist << " which is less than the"
776 << std::endl;
777 oomph_info << "permitted distance of " << epsilon << std::endl
778 << std::endl;
779 oomph_info << "The offending nodes are located at: " << std::endl;
780 for (unsigned i = 0; i < dim; i++)
781 {
782 oomph_info << nod1_pt->x(i) << " ";
783 }
784 if (nod1_pt->is_a_copy() || nod2_pt->is_a_copy())
785 {
786 oomph_info << "\n\n[NOTE: message issued as diagonistic rather "
787 "than an error\n"
788 << " because at least one of the nodes is a copy; you "
789 "may still\n"
790 << " want to check this out. BACKGROUND: Copied nodes "
791 "share the same Data but\n"
792 << " will, in general, have different spatial "
793 "positions (e.g. when used\n"
794 << " as periodic nodes); however there are cases when "
795 "they are located\n"
796 << " at the same spatial position (e.g. in "
797 "oomph-lib's annular mesh which\n"
798 << " is a rolled-around version of the rectangular "
799 "quadmesh). In such cases,\n"
800 << " the nodes could have been deleted and completely "
801 "replaced by \n"
802 << " pointers to existing nodes, but may have been "
803 "left there for convenience\n"
804 << " or out of laziness...]\n";
805 }
806 else
807 {
808 failed = true;
809 }
810 oomph_info << std::endl << std::endl;
811 }
812 }
813 }
814 if (failed) return 1;
815
816 // If we made it to here, we must have passed the test.
817 oomph_info << "...done: Test passed!" << std::endl << std::endl;
818 return 0;
819 }
820
821 /// Prune nodes. Nodes that have been marked as obsolete are removed
822 /// from the mesh (and its boundary-node scheme). Returns vector
823 /// of pointers to deleted nodes.
825
826 /// Return number of boundaries
827 unsigned nboundary() const
828 {
829 return Boundary_node_pt.size();
830 }
831
832 /// Return number of nodes on a particular boundary
833 unsigned long nboundary_node(const unsigned& ibound) const
834 {
835 return Boundary_node_pt[ibound].size();
836 }
837
838
839 /// Return pointer to e-th finite element on boundary b
841 const unsigned& e) const
842 {
843#ifdef PARANOID
845 {
846 throw OomphLibError("Lookup scheme for elements next to boundary "
847 "hasn't been set up yet!\n",
848 OOMPH_CURRENT_FUNCTION,
849 OOMPH_EXCEPTION_LOCATION);
850 }
851#endif
852 return Boundary_element_pt[b][e];
853 }
854
855
856 /// Find a node not on any boundary in mesh_pt (useful for pinning
857 /// a single node in a purely Neumann problem so that it is fully
858 /// determined).
860 {
861 for (unsigned nd = 0, nnd = nnode(); nd < nnd; nd++)
862 {
863 if (!(node_pt(nd)->is_on_boundary()))
864 {
865 return node_pt(nd);
866 }
867 }
868
869 std::ostringstream error_msg;
870 error_msg << "No non-boundary nodes in the mesh.";
871 throw OomphLibError(
872 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
873 // Never get here!
874 return 0;
875 }
876
877 /// Return number of finite elements that are adjacent to boundary b
878 unsigned nboundary_element(const unsigned& b) const
879 {
880#ifdef PARANOID
882 {
883 throw OomphLibError("Lookup scheme for elements next to boundary "
884 "hasn't been set up yet!\n",
885 OOMPH_CURRENT_FUNCTION,
886 OOMPH_EXCEPTION_LOCATION);
887 }
888#endif
889 return Boundary_element_pt[b].size();
890 }
891
892
893 /// For the e-th finite element on boundary b, return int to indicate
894 /// the face_index of the face adjacent to the boundary. This is consistent
895 /// with input required during the generation of FaceElements.
896 int face_index_at_boundary(const unsigned& b, const unsigned& e) const
897 {
898#ifdef PARANOID
900 {
901 throw OomphLibError("Lookup scheme for elements next to boundary "
902 "hasn't been set up yet!\n",
903 OOMPH_CURRENT_FUNCTION,
904 OOMPH_EXCEPTION_LOCATION);
905 }
906#endif
907 return Face_index_at_boundary[b][e];
908 }
909
910 /// Dump the data in the mesh into a file for restart
911 virtual void dump(std::ofstream& dump_file,
912 const bool& use_old_ordering = true) const;
913
914 /// Dump the data in the mesh into a file for restart
915 void dump(const std::string& dump_file_name,
916 const bool& use_old_ordering = true) const
917 {
918 std::ofstream dump_stream(dump_file_name.c_str());
919#ifdef PARANOID
920 if (!dump_stream.is_open())
921 {
922 std::string err = "Couldn't open file " + dump_file_name;
923 throw OomphLibError(
924 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
925 }
926#endif
927 dump(dump_stream, use_old_ordering);
928 }
929
930 /// Read solution from restart file
931 virtual void read(std::ifstream& restart_file);
932
933
934 /// Output in paraview format into specified file. Breaks up each
935 /// element into sub-elements for plotting purposes. We assume
936 /// that all elements are of the same type (fct will break
937 /// break (in paranoid mode) if paraview output fcts of the
938 /// elements are inconsistent).
939 void output_paraview(std::ofstream& file_out, const unsigned& nplot) const;
940
941 /// Output in paraview format into specified file. Breaks up each
942 /// element into sub-elements for plotting purposes. We assume
943 /// that all elements are of the same type (fct will break
944 /// break (in paranoid mode) if paraview output fcts of the
945 /// elements are inconsistent).
947 std::ofstream& file_out,
948 const unsigned& nplot,
949 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const;
950
951 /// Output in paraview format into specified file. Breaks up each
952 /// element into sub-elements for plotting purposes. We assume
953 /// that all elements are of the same type (fct will break
954 /// break (in paranoid mode) if paraview output fcts of the
955 /// elements are inconsistent).
957 std::ofstream& file_out,
958 const unsigned& nplot,
959 const double& time,
961
962 /// Output for all elements
963 void output(std::ostream& outfile);
964
965 /// Output at f(n_plot) points in each element
966 void output(std::ostream& outfile, const unsigned& n_plot);
967
968 /// Output for all elements (C-style output)
969 void output(FILE* file_pt);
970
971 /// Output at f(n_plot) points in each element (C-style output)
972 void output(FILE* file_pt, const unsigned& nplot);
973
974 /// Output for all elements
975 void output(const std::string& output_filename)
976 {
977 std::ofstream outfile;
978 outfile.open(output_filename.c_str());
979 output(outfile);
980 outfile.close();
981 }
982
983 /// Output at f(n_plot) points in each element
984 void output(const std::string& output_filename, const unsigned& n_plot)
985 {
986 std::ofstream outfile;
987 outfile.open(output_filename.c_str());
988 output(outfile, n_plot);
989 outfile.close();
990 }
991
992 /// Output a given Vector function at f(n_plot) points in each element
993 void output_fct(std::ostream& outfile,
994 const unsigned& n_plot,
996
997 /// Output a given time-dep. Vector function at f(n_plot) points in
998 /// each element
999 void output_fct(std::ostream& outfile,
1000 const unsigned& n_plot,
1001 const double& time,
1003
1004 /// Output the nodes on the boundaries (into separate tecplot zones)
1005 void output_boundaries(std::ostream& outfile);
1006
1007 /// Output the nodes on the boundaries (into separate tecplot zones).
1008 /// Specify filename
1009 void output_boundaries(const std::string& output_filename)
1010 {
1011 std::ofstream outfile;
1012 outfile.open(output_filename.c_str());
1013 output_boundaries(outfile);
1014 outfile.close();
1015 }
1016
1017 /// Assign initial values for an impulsive start
1019
1020 /// Shift time-dependent data along for next timestep:
1021 /// Deal with nodal Data/positions and the element's internal
1022 /// Data
1023 void shift_time_values();
1024
1025
1026 /// Calculate predictions for all Data and positions associated
1027 /// with the mesh, usually used in adaptive time-stepping.
1028 void calculate_predictions();
1029
1030 /// Set the timestepper associated with all nodal and elemental
1031 /// data stored in the mesh.
1033 TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data)
1034 {
1035 this->set_nodal_time_stepper(time_stepper_pt, preserve_existing_data);
1036 this->set_elemental_internal_time_stepper(time_stepper_pt,
1037 preserve_existing_data);
1038 }
1039
1040 /// Function that can be used to set any additional timestepper data
1041 /// stored at the Mesh (as opposed to nodal and elemental) levels. This
1042 /// is virtual so that it can be overloaded in the appropriate Meshes.
1043 /// Examples include the SpineMeshes and adaptive triangle and tet meshes
1044 virtual void set_mesh_level_time_stepper(
1045 TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
1046
1047 /// Set consistent values for pinned data in continuation
1049 ContinuationStorageScheme* const& continuation_stepper_pt);
1050
1051 /// Does the double pointer correspond to any mesh data
1052 bool does_pointer_correspond_to_mesh_data(double* const& parameter_pt);
1053
1054 /// Set the timestepper associated with the nodal data in the mesh
1055 void set_nodal_time_stepper(TimeStepper* const& time_stepper_pt,
1056 const bool& preserve_existing_data);
1057
1058 /// Set the timestepper associated with the internal data stored
1059 /// within elements in the meah
1061 TimeStepper* const& time_stepper_pt, const bool& preserve_existing_data);
1062
1063 /// Compute norm of solution by summing contributions of
1064 /// compute_norm(...) for all constituent elements in the mesh.
1065 /// What that norm means depends on what's defined in the element's
1066 /// function; may need to take the square root afterwards if the elements
1067 /// compute the square of the L2 norm, say.
1068 virtual void compute_norm(double& norm)
1069 {
1070 // Initialse the norm
1071 norm = 0.0;
1072
1073 // Per-element norm
1074 double el_norm = 0;
1075
1076 // Loop over the elements
1077 unsigned long n_element = Element_pt.size();
1078 for (unsigned long e = 0; e < n_element; e++)
1079 {
1081
1082#ifdef OOMPH_HAS_MPI
1083 // Compute error for each non-halo element
1084 if (!(el_pt->is_halo()))
1085#endif
1086 {
1087 el_pt->compute_norm(el_norm);
1088 }
1089 norm += el_norm;
1090 }
1091 }
1092
1093
1094 /// Compute norm of solution by summing contributions of
1095 /// compute_norm(...) for all constituent elements in the mesh.
1096 /// What that norm means depends on what's defined in the element's
1097 /// function; may need to take the square root afterwards if the elements
1098 /// compute the square of the L2 norm, say.
1099 virtual void compute_norm(Vector<double>& norm)
1100 {
1101 // How many unknowns are there?
1102 unsigned n_entry = norm.size();
1103
1104 // Initialse the norm
1105 norm.initialise(0.0);
1106
1107 // Per-element norm
1108 Vector<double> el_norm(n_entry, 0.0);
1109
1110 // How many elements are there?
1111 unsigned long n_element = Element_pt.size();
1112
1113 // Loop over the elements
1114 for (unsigned long e = 0; e < n_element; e++)
1115 {
1116 // Get a pointer to the e-th generalised element
1118
1119#ifdef OOMPH_HAS_MPI
1120 // Compute error for each non-halo element
1121 if (!(el_pt->is_halo()))
1122#endif
1123 {
1124 // Compute the elemental norm
1125 el_pt->compute_norm(el_norm);
1126 }
1127
1128 // Loop over the norm vector entries
1129 for (unsigned i = 0; i < n_entry; i++)
1130 {
1131 // Update the norm of the i-th component
1132 norm[i] += el_norm[i];
1133 }
1134 } // for (unsigned long e=0;e<n_element;e++)
1135 } // End of compute_norm
1136
1137
1138 /// Plot error when compared against a given exact solution.
1139 /// Also returns the norm of the error and that of the exact solution
1140 virtual void compute_error(
1141 std::ostream& outfile,
1143 const double& time,
1144 double& error,
1145 double& norm)
1146 {
1147 // Initialse the norm and error
1148 norm = 0.0;
1149 error = 0.0;
1150 // Per-element norm and error
1151 double el_error, el_norm;
1152
1153 // Loop over the elements
1154 unsigned long Element_pt_range = Element_pt.size();
1155 for (unsigned long e = 0; e < Element_pt_range; e++)
1156 {
1157 // Try to cast to FiniteElement
1158 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1159 if (el_pt == 0)
1160 {
1161 throw OomphLibError(
1162 "Can't execute compute_error(...) for non FiniteElements",
1163 OOMPH_CURRENT_FUNCTION,
1164 OOMPH_EXCEPTION_LOCATION);
1165 }
1166
1167 // Reset elemental errors and norms
1168 el_norm = 0.0;
1169 el_error = 0.0;
1170#ifdef OOMPH_HAS_MPI
1171 // Compute error for each non-halo element
1172 if (!(el_pt->is_halo()))
1173#endif
1174 {
1175 el_pt->compute_error(outfile, exact_soln_pt, time, el_error, el_norm);
1176 }
1177 // Add each element error to the global errors
1178 norm += el_norm;
1179 error += el_error;
1180 }
1181 }
1182
1183 /// Plot error when compared against a given time-depdendent
1184 /// exact solution. Also returns the norm of the error and
1185 /// that of the exact solution
1186 virtual void compute_error(
1187 std::ostream& outfile,
1189 double& error,
1190 double& norm)
1191 {
1192 // Initialise norm and error
1193 norm = 0.0;
1194 error = 0.0;
1195 // Per-element norm and error
1196 double el_error, el_norm;
1197
1198 // Loop over the elements
1199 unsigned long Element_pt_range = Element_pt.size();
1200 for (unsigned long e = 0; e < Element_pt_range; e++)
1201 {
1202 // Try to cast to FiniteElement
1203 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1204 if (el_pt == 0)
1205 {
1206 throw OomphLibError(
1207 "Can't execute compute_error(...) for non FiniteElements",
1208 OOMPH_CURRENT_FUNCTION,
1209 OOMPH_EXCEPTION_LOCATION);
1210 }
1211 // Reset elemental errors and norms
1212 el_norm = 0.0;
1213 el_error = 0.0;
1214 // Calculate the elemental errors for each non-halo element
1215#ifdef OOMPH_HAS_MPI
1216 if (!(el_pt->is_halo()))
1217#endif
1218 {
1219 el_pt->compute_error(outfile, exact_soln_pt, el_error, el_norm);
1220 }
1221 // Add each elemental error to the global error
1222 norm += el_norm;
1223 error += el_error;
1224 }
1225 }
1226
1227
1228 /// Plot error when compared against a given time-dependent
1229 /// exact solution. Also returns the norm of the error and
1230 /// that of the exact solution
1231 virtual void compute_error(
1233 double& error,
1234 double& norm)
1235 {
1236 // Initialise norm and error
1237 norm = 0.0;
1238 error = 0.0;
1239
1240 // Per-element norm and error
1241 double el_error, el_norm;
1242
1243 // Loop over the elements
1244 unsigned long Element_pt_range = Element_pt.size();
1245 for (unsigned long e = 0; e < Element_pt_range; e++)
1246 {
1247 // Try to cast to FiniteElement
1248 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1249 if (el_pt == 0)
1250 {
1251 throw OomphLibError(
1252 "Can't execute compute_error(...) for non FiniteElements",
1253 OOMPH_CURRENT_FUNCTION,
1254 OOMPH_EXCEPTION_LOCATION);
1255 }
1256
1257 // Reset elemental errors and norms
1258 el_norm = 0.0;
1259 el_error = 0.0;
1260
1261 // Calculate the elemental errors for each non-halo element
1262#ifdef OOMPH_HAS_MPI
1263 if (!(el_pt->is_halo()))
1264#endif
1265 {
1266 el_pt->compute_error(exact_soln_pt, el_error, el_norm);
1267 }
1268
1269 // Add each elemental error to the global error
1270 norm += el_norm;
1271 error += el_error;
1272 }
1273 } // End of compute_error
1274
1275
1276 /// Plot error when compared against a given time-dependent
1277 /// exact solution. Also returns the norm of the error and
1278 /// that of the exact solution
1279 virtual void compute_error(
1281 Vector<double>& error,
1282 Vector<double>& norm)
1283 {
1284 // Initialise norm vector entries
1285 norm.initialise(0.0);
1286
1287 // Initialise error vector entries
1288 error.initialise(0.0);
1289
1290 // Norm vector size
1291 unsigned n_norm = norm.size();
1292
1293 // Error vector size
1294 unsigned n_error = error.size();
1295
1296 // Per-element norm and error
1297 Vector<double> el_norm(n_norm, 0.0);
1298
1299 // Per-element norm and error
1300 Vector<double> el_error(n_error, 0.0);
1301
1302 // How many elements are there?
1303 unsigned long element_pt_range = Element_pt.size();
1304
1305 // Loop over the elements
1306 for (unsigned long e = 0; e < element_pt_range; e++)
1307 {
1308 // Try to cast to FiniteElement
1309 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1310 if (el_pt == 0)
1311 {
1312 throw OomphLibError(
1313 "Can't execute compute_error(...) for non FiniteElements",
1314 OOMPH_CURRENT_FUNCTION,
1315 OOMPH_EXCEPTION_LOCATION);
1316 }
1317
1318 // Re-initialise norm vector entries
1319 el_norm.initialise(0.0);
1320
1321 // Re-initialise error vector entries
1322 el_error.initialise(0.0);
1323
1324 // Calculate the elemental errors for each non-halo element
1325#ifdef OOMPH_HAS_MPI
1326 if (!(el_pt->is_halo()))
1327#endif
1328 {
1329 el_pt->compute_error(exact_soln_pt, el_error, el_norm);
1330 }
1331
1332 // Add each elemental norm contribution to the global norm
1333 for (unsigned i = 0; i < n_norm; i++)
1334 {
1335 norm[i] += el_norm[i];
1336 }
1337
1338 // Add each elemental error contribution to the global error
1339 for (unsigned i = 0; i < n_error; i++)
1340 {
1341 error[i] += el_error[i];
1342 }
1343 } // for (unsigned long e=0;e<element_pt_range;e++)
1344 } // End of compute_error
1345
1346
1347 /// Plot error when compared against a given time-depdendent
1348 /// exact solution. Also returns the norm of the error and
1349 /// that of the exact solution. Version with vectors of norms and errors so
1350 /// that different variables' norms and errors can be returned individually
1351 virtual void compute_error(
1352 std::ostream& outfile,
1354 const double& time,
1355 Vector<double>& error,
1356 Vector<double>& norm)
1357 {
1358 // Initialise norm and error
1359 unsigned n_error = error.size();
1360 unsigned n_norm = norm.size();
1361 for (unsigned i = 0; i < n_error; i++)
1362 {
1363 error[i] = 0.0;
1364 }
1365 for (unsigned i = 0; i < n_norm; i++)
1366 {
1367 norm[i] = 0.0;
1368 }
1369 // Per-element norm and error
1370 Vector<double> el_error(n_error), el_norm(n_norm);
1371
1372 // Loop over the elements
1373 unsigned long Element_pt_range = Element_pt.size();
1374 for (unsigned long e = 0; e < Element_pt_range; e++)
1375 {
1376 // Try to cast to FiniteElement
1377 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1378 if (el_pt == 0)
1379 {
1380 throw OomphLibError(
1381 "Can't execute compute_error(...) for non FiniteElements",
1382 OOMPH_CURRENT_FUNCTION,
1383 OOMPH_EXCEPTION_LOCATION);
1384 }
1385 // Reset elemental errors and norms
1386 for (unsigned i = 0; i < n_error; i++)
1387 {
1388 el_error[i] = 0.0;
1389 }
1390 for (unsigned i = 0; i < n_norm; i++)
1391 {
1392 el_norm[i] = 0.0;
1393 }
1394 // Calculate the elemental errors for each non-halo element
1395#ifdef OOMPH_HAS_MPI
1396 if (!(el_pt->is_halo()))
1397#endif
1398 {
1399 el_pt->compute_error(outfile, exact_soln_pt, time, el_error, el_norm);
1400 }
1401 // Add each elemental error to the global error
1402 for (unsigned i = 0; i < n_error; i++)
1403 {
1404 error[i] += el_error[i];
1405 }
1406 for (unsigned i = 0; i < n_norm; i++)
1407 {
1408 norm[i] += el_norm[i];
1409 }
1410 }
1411 }
1412
1413 /// Plot error when compared against a given time-depdendent
1414 /// exact solution. Also returns the norm of the error and
1415 /// that of the exact solution. Version with vectors of norms and errors so
1416 /// that different variables' norms and errors can be returned individually
1417 virtual void compute_error(
1418 std::ostream& outfile,
1420 Vector<double>& error,
1421 Vector<double>& norm)
1422 {
1423 // Initialise norm and error
1424 unsigned n_error = error.size();
1425 unsigned n_norm = norm.size();
1426 for (unsigned i = 0; i < n_error; i++)
1427 {
1428 error[i] = 0.0;
1429 }
1430 for (unsigned i = 0; i < n_norm; i++)
1431 {
1432 norm[i] = 0.0;
1433 }
1434 // Per-element norm and error
1435 Vector<double> el_error(n_error), el_norm(n_norm);
1436
1437 // Loop over the elements
1438 unsigned long Element_pt_range = Element_pt.size();
1439 for (unsigned long e = 0; e < Element_pt_range; e++)
1440 {
1441 // Try to cast to FiniteElement
1442 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1443 if (el_pt == 0)
1444 {
1445 throw OomphLibError(
1446 "Can't execute compute_error(...) for non FiniteElements",
1447 OOMPH_CURRENT_FUNCTION,
1448 OOMPH_EXCEPTION_LOCATION);
1449 }
1450 // Reset elemental errors and norms
1451 for (unsigned i = 0; i < n_error; i++)
1452 {
1453 el_error[i] = 0.0;
1454 }
1455 for (unsigned i = 0; i < n_norm; i++)
1456 {
1457 el_norm[i] = 0.0;
1458 }
1459 // Calculate the elemental errors for each non-halo element
1460#ifdef OOMPH_HAS_MPI
1461 if (!(el_pt->is_halo()))
1462#endif
1463 {
1464 el_pt->compute_error(outfile, exact_soln_pt, el_error, el_norm);
1465 }
1466 // Add each elemental error to the global error
1467 for (unsigned i = 0; i < n_error; i++)
1468 {
1469 error[i] += el_error[i];
1470 }
1471 for (unsigned i = 0; i < n_norm; i++)
1472 {
1473 norm[i] += el_norm[i];
1474 }
1475 }
1476 }
1477
1478 /// Returns the norm of the error and that of the exact solution
1479 virtual void compute_error(
1481 const double& time,
1482 double& error,
1483 double& norm)
1484 {
1485 // Initialse the norm and error
1486 norm = 0.0;
1487 error = 0.0;
1488 // Per-element norm and error
1489 double el_error, el_norm;
1490
1491 // Loop over the elements
1492 unsigned long Element_pt_range = Element_pt.size();
1493 for (unsigned long e = 0; e < Element_pt_range; e++)
1494 {
1495 // Try to cast to FiniteElement
1496 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1497 if (el_pt == 0)
1498 {
1499 throw OomphLibError(
1500 "Can't execute compute_error(...) for non FiniteElements",
1501 OOMPH_CURRENT_FUNCTION,
1502 OOMPH_EXCEPTION_LOCATION);
1503 }
1504
1505 // Reset elemental errors and norms
1506 el_norm = 0.0;
1507 el_error = 0.0;
1508#ifdef OOMPH_HAS_MPI
1509 // Compute error for each non-halo element
1510 if (!(el_pt->is_halo()))
1511#endif
1512 {
1513 el_pt->compute_error(exact_soln_pt, time, el_error, el_norm);
1514 }
1515 // Add each element error to the global errors
1516 norm += el_norm;
1517 error += el_error;
1518 }
1519 }
1520
1521
1522 /// Returns the norm of the error and that of the exact solution.
1523 /// Version with vectors of norms and errors so that different variables'
1524 /// norms and errors can be returned individually
1525 virtual void compute_error(
1527 const double& time,
1528 Vector<double>& error,
1529 Vector<double>& norm)
1530 {
1531 // Initialise norm and error
1532 unsigned n_error = error.size();
1533 unsigned n_norm = norm.size();
1534 for (unsigned i = 0; i < n_error; i++)
1535 {
1536 error[i] = 0.0;
1537 }
1538 for (unsigned i = 0; i < n_norm; i++)
1539 {
1540 norm[i] = 0.0;
1541 }
1542 // Per-element norm and error
1543 Vector<double> el_error(n_error), el_norm(n_norm);
1544
1545 // Loop over the elements
1546 unsigned long Element_pt_range = Element_pt.size();
1547 for (unsigned long e = 0; e < Element_pt_range; e++)
1548 {
1549 // Try to cast to FiniteElement
1550 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
1551 if (el_pt == 0)
1552 {
1553 throw OomphLibError(
1554 "Can't execute compute_error(...) for non FiniteElements",
1555 OOMPH_CURRENT_FUNCTION,
1556 OOMPH_EXCEPTION_LOCATION);
1557 }
1558 // Reset elemental errors and norms
1559 for (unsigned i = 0; i < n_error; i++)
1560 {
1561 el_error[i] = 0.0;
1562 }
1563 for (unsigned i = 0; i < n_norm; i++)
1564 {
1565 el_norm[i] = 0.0;
1566 }
1567 // Calculate the elemental errors for each non-halo element
1568#ifdef OOMPH_HAS_MPI
1569 if (!(el_pt->is_halo()))
1570#endif
1571 {
1572 el_pt->compute_error(exact_soln_pt, time, el_error, el_norm);
1573 }
1574 // Add each elemental error to the global error
1575 for (unsigned i = 0; i < n_error; i++)
1576 {
1577 error[i] += el_error[i];
1578 }
1579 for (unsigned i = 0; i < n_norm; i++)
1580 {
1581 norm[i] += el_norm[i];
1582 }
1583 } // for (unsigned long e=0;e<Element_pt_range;e++)
1584 } // End of compute_error
1585
1586
1587 /// Boolean to indicate if Mesh has been distributed
1589 {
1590#ifdef OOMPH_HAS_MPI
1591 return communicator_pt() != 0;
1592#else
1593 // Can't be distributed without mpi
1594 return false;
1595#endif
1596 }
1597
1598 /// Read-only access fct to communicator (Null if mesh is not distributed,
1599 /// i.e. if we don't have mpi).
1601 {
1602#ifdef OOMPH_HAS_MPI
1603 return Comm_pt;
1604#else
1605 return 0;
1606#endif
1607 }
1608
1609
1610#ifdef OOMPH_HAS_MPI
1611
1612 /// Function to set communicator (mesh is assumed to be distributed if the
1613 /// communicator pointer is non-null). Only defined if mpi is enabled
1614 /// becaus Comm_pt itself is only defined when mpi is enabled.
1616 {
1617 Comm_pt = comm_pt;
1618 }
1619
1620 /// Call this function to keep all the elements as halo elements
1622 {
1624 }
1625
1626 /// Calll this function to unset the flag that keeps all elements
1627 /// in the mesh as halo elements
1629 {
1631 }
1632
1633 /// Distribute the problem and doc; make this virtual to allow
1634 /// overloading for particular meshes where further work is required.
1635 /// Add to vector of pointers to deleted elements.
1636 virtual void distribute(OomphCommunicator* comm_pt,
1637 const Vector<unsigned>& element_domain,
1638 Vector<GeneralisedElement*>& deleted_element_pt,
1639 DocInfo& doc_info,
1640 const bool& report_stats,
1641 const bool& overrule_keep_as_halo_element_status);
1642
1643 /// Distribute the problem
1644 /// Add to vector of pointers to deleted elements.
1646 const Vector<unsigned>& element_domain,
1647 Vector<GeneralisedElement*>& deleted_element_pt,
1648 const bool& report_stats = false)
1649 {
1650 DocInfo doc_info;
1651 doc_info.disable_doc();
1652 bool overrule_keep_as_halo_element_status = false;
1653 return distribute(comm_pt,
1654 element_domain,
1655 deleted_element_pt,
1656 doc_info,
1657 report_stats,
1658 overrule_keep_as_halo_element_status);
1659 }
1660
1661 /// (Irreversibly) prune halo(ed) elements and nodes, usually
1662 /// after another round of refinement, to get rid of
1663 /// excessively wide halo layers. Note that the current
1664 /// mesh will be now regarded as the base mesh and no unrefinement
1665 /// relative to it will be possible once this function
1666 /// has been called.
1668 Vector<GeneralisedElement*>& deleted_element_pt,
1669 const bool& report_stats = false)
1670 {
1671 DocInfo doc_info;
1672 doc_info.disable_doc();
1673 prune_halo_elements_and_nodes(deleted_element_pt, doc_info, report_stats);
1674 }
1675
1676
1677 /// (Irreversibly) prune halo(ed) elements and nodes, usually
1678 /// after another round of refinement, to get rid of
1679 /// excessively wide halo layers. Note that the current
1680 /// mesh will be now regarded as the base mesh and no unrefinement
1681 /// relative to it will be possible once this function
1682 /// has been called.
1684 Vector<GeneralisedElement*>& deleted_element_pt,
1685 DocInfo& doc_info,
1686 const bool& report_stats);
1687
1688 /// Get efficiency of mesh distribution: In an ideal distribution
1689 /// without halo overhead, each processor would only hold its own
1690 /// elements. Efficieny per processor = (number of non-halo elements)/
1691 /// (total number of elements).
1692 void get_efficiency_of_mesh_distribution(double& av_efficiency,
1693 double& max_efficiency,
1694 double& min_efficiency);
1695
1696 /// Doc the mesh distribution, to be processed with tecplot macros
1697 void doc_mesh_distribution(DocInfo& doc_info);
1698
1699 /// Check halo and shared schemes on the mesh
1700 void check_halo_schemes(DocInfo& doc_info,
1701 double& max_permitted_error_for_halo_check);
1702
1703 /// Classify the halo and haloed nodes in the mesh. Virtual
1704 /// so it can be overloaded to perform additional functionality
1705 /// (such as synchronising hanging nodes) in refineable meshes, say.
1706 virtual void classify_halo_and_haloed_nodes(DocInfo& doc_info,
1707 const bool& report_stats);
1708
1709 /// Classify the halo and haloed nodes in the mesh. Virtual
1710 /// so it can be overloaded to perform additional functionality
1711 /// (such as synchronising hanging nodes) in refineable meshes, say.
1713 const bool& report_stats = false)
1714 {
1715 DocInfo doc_info;
1716 doc_info.disable_doc();
1717 classify_halo_and_haloed_nodes(doc_info, report_stats);
1718 }
1719
1720 /// Synchronise shared node lookup schemes to cater for the
1721 /// the case where:
1722 /// (1) a certain node on the current processor is halo with proc p
1723 /// (i.e. its non-halo counterpart lives on processor p)
1724 /// (2) that node is also exists (also as a halo) on another processor
1725 /// (q, say) where its non-halo counter part is also known to be
1726 /// on processor p.
1727 /// However, without calling this function the current processor does not
1728 /// necessarily know that it shares a node with processor q. This
1729 /// information can be required, e.g. when synchronising hanging node
1730 /// schemes over all processors.
1731 void synchronise_shared_nodes(const bool& report_stats);
1732
1733
1734 /// Get all the halo data stored in the mesh and add pointers to
1735 /// the data to the map, indexed by global equation number
1736 void get_all_halo_data(std::map<unsigned, double*>& map_of_halo_data);
1737
1738 /// Return vector of halo elements in this Mesh
1739 /// whose non-halo counterpart is held on processor p.
1741 {
1742 // Prepare vector
1744
1745 // Loop over all root halo elements
1746 unsigned nelem = nroot_halo_element(p);
1747 for (unsigned e = 0; e < nelem; e++)
1748 {
1750
1751 // Is it a refineable element?
1752 RefineableElement* ref_el_pt = dynamic_cast<RefineableElement*>(el_pt);
1753 if (ref_el_pt != 0)
1754 {
1755 // Vector of pointers to leaves in tree emanating from
1756 // current root halo element
1757 Vector<Tree*> leaf_pt;
1758 ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
1759
1760 // Loop over leaves and add their objects (the finite elements)
1761 // to vector
1762 unsigned nleaf = leaf_pt.size();
1763 for (unsigned l = 0; l < nleaf; l++)
1764 {
1765 vec_el_pt.push_back(leaf_pt[l]->object_pt());
1766 }
1767 }
1768 else
1769 {
1770 vec_el_pt.push_back(el_pt);
1771 }
1772 }
1773 return vec_el_pt;
1774 }
1775
1776
1777 /// Return vector of haloed elements in this Mesh
1778 /// whose haloing counterpart is held on processor p.
1780 {
1781 // Prepare vector
1783
1784 // Loop over all root haloed elements
1785 unsigned nelem = nroot_haloed_element(p);
1786 for (unsigned e = 0; e < nelem; e++)
1787 {
1789
1790 // Is it a refineable element?
1791 RefineableElement* ref_el_pt = dynamic_cast<RefineableElement*>(el_pt);
1792 if (ref_el_pt != 0)
1793 {
1794 // Vector of pointers to leaves in tree emanating from
1795 // current root haloed element
1796 Vector<Tree*> leaf_pt;
1797 ref_el_pt->tree_pt()->stick_leaves_into_vector(leaf_pt);
1798
1799 // Loop over leaves and add their objects (the finite elements)
1800 // to vector
1801 unsigned nleaf = leaf_pt.size();
1802 for (unsigned l = 0; l < nleaf; l++)
1803 {
1804 vec_el_pt.push_back(leaf_pt[l]->object_pt());
1805 }
1806 }
1807 else
1808 {
1809 vec_el_pt.push_back(el_pt);
1810 }
1811 }
1812 return vec_el_pt;
1813 }
1814
1815 /// Total number of non-halo elements in this mesh (Costly call
1816 /// computes result on the fly)
1818 {
1819 unsigned count = 0;
1820 unsigned n = nelement();
1821 for (unsigned e = 0; e < n; e++)
1822 {
1823 if (!(element_pt(e)->is_halo())) count++;
1824 }
1825 return count;
1826 }
1827
1828
1829 /// Total number of root halo elements in this Mesh
1831 {
1832 unsigned n = 0;
1833 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
1834 Root_halo_element_pt.begin();
1835 it != Root_halo_element_pt.end();
1836 it++)
1837 {
1838 n += it->second.size();
1839 }
1840 return n;
1841 }
1842
1843
1844 /// Number of root halo elements in this Mesh whose non-halo
1845 /// counterpart is held on processor p.
1846 unsigned nroot_halo_element(const unsigned& p)
1847 {
1848 return Root_halo_element_pt[p].size();
1849 }
1850
1851
1852 /// Vector of pointers to root halo elements in this Mesh
1853 /// whose non-halo counterpart is held on processor p.
1855 {
1856 return Root_halo_element_pt[p];
1857 }
1858
1859
1860 /// Access fct to the e-th root halo element in this Mesh
1861 /// whose non-halo counterpart is held on processor p.
1863 const unsigned& e)
1864 {
1865 return Root_halo_element_pt[p][e];
1866 }
1867
1868
1869 /// Add root halo element whose non-halo counterpart is held
1870 /// on processor p to this Mesh.
1871 void add_root_halo_element_pt(const unsigned& p, GeneralisedElement*& el_pt)
1872 {
1873 Root_halo_element_pt[p].push_back(el_pt);
1874 el_pt->set_halo(p);
1875 }
1876
1877 /// Total number of halo nodes in this Mesh
1878 unsigned nhalo_node()
1879 {
1880 unsigned n = 0;
1881 for (std::map<unsigned, Vector<Node*>>::iterator it =
1882 Halo_node_pt.begin();
1883 it != Halo_node_pt.end();
1884 it++)
1885 {
1886 n += it->second.size();
1887 }
1888 return n;
1889 }
1890
1891 /// Number of halo nodes in this Mesh whose non-halo counterpart
1892 /// is held on processor p.
1893 unsigned nhalo_node(const unsigned& p)
1894 {
1895 // Memory saving version of: return Halo_node_pt[p].size();
1896 std::map<unsigned, Vector<Node*>>::iterator it = Halo_node_pt.find(p);
1897 if (it == Halo_node_pt.end())
1898 {
1899 return 0;
1900 }
1901 return (*it).second.size();
1902 }
1903
1904 /// Add halo node whose non-halo counterpart is held
1905 /// on processor p to the storage scheme for halo nodes.
1906 void add_halo_node_pt(const unsigned& p, Node*& nod_pt)
1907 {
1908 Halo_node_pt[p].push_back(nod_pt);
1909 }
1910
1911
1912 /// Access fct to the j-th halo node in this Mesh
1913 /// whose non-halo counterpart is held on processor p.
1914 Node* halo_node_pt(const unsigned& p, const unsigned& j)
1915 {
1916 return Halo_node_pt[p][j];
1917 }
1918
1919
1920 /// Total number of root haloed elements in this Mesh
1922 {
1923 unsigned n = 0;
1924 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
1925 Root_haloed_element_pt.begin();
1926 it != Root_haloed_element_pt.end();
1927 it++)
1928 {
1929 n += it->second.size();
1930 }
1931 return n;
1932 }
1933
1934 /// Number of root haloed elements in this Mesh whose non-halo
1935 /// counterpart is held on processor p.
1936 unsigned nroot_haloed_element(const unsigned& p)
1937 {
1938 // Memory saving version of: return Root_haloed_element_pt[p].size();
1939 std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
1940 Root_haloed_element_pt.find(p);
1941 if (it == Root_haloed_element_pt.end())
1942 {
1943 return 0;
1944 }
1945 return (*it).second.size();
1946 }
1947
1948
1949 /// Vector of pointers to root haloed elements in this Mesh
1950 /// whose non-halo counterpart is held on processor p.
1952 {
1953 // Memory saving version of: return Root_haloed_element_pt[p];
1954 std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
1955 Root_haloed_element_pt.find(p);
1956 if (it == Root_haloed_element_pt.end())
1957 {
1959 return tmp;
1960 }
1961 return (*it).second;
1962 }
1963
1964 /// Access fct to the e-th root haloed element in this Mesh
1965 /// whose non-halo counterpart is held on processor p.
1967 const unsigned& e)
1968 {
1969 return Root_haloed_element_pt[p][e];
1970 }
1971
1972 /// Add root haloed element whose non-halo counterpart is held
1973 /// on processor p to the storage scheme for haloed elements.
1974 /// Note: This does not add the element to the storage scheme
1975 /// for elements as it's understood to naturally live on this
1976 /// processor anyway!
1977 void add_root_haloed_element_pt(const unsigned& p,
1978 GeneralisedElement*& el_pt)
1979 {
1980 Root_haloed_element_pt[p].push_back(el_pt);
1981 }
1982
1983
1984 /// Total number of haloed nodes in this Mesh
1985 unsigned nhaloed_node()
1986 {
1987 unsigned n = 0;
1988 for (std::map<unsigned, Vector<Node*>>::iterator it =
1989 Haloed_node_pt.begin();
1990 it != Haloed_node_pt.end();
1991 it++)
1992 {
1993 n += it->second.size();
1994 }
1995 return n;
1996 }
1997
1998
1999 /// Number of haloed nodes in this Mesh whose haloed counterpart
2000 /// is held on processor p.
2001 unsigned nhaloed_node(const unsigned& p)
2002 {
2003 // Memory saving version of: return Haloed_node_pt[p].size();
2004 std::map<unsigned, Vector<Node*>>::iterator it = Haloed_node_pt.find(p);
2005 if (it == Haloed_node_pt.end())
2006 {
2007 return 0;
2008 }
2009 return (*it).second.size();
2010 }
2011
2012 /// Access fct to the j-th haloed node in this Mesh
2013 /// whose halo counterpart is held on processor p.
2014 Node* haloed_node_pt(const unsigned& p, const unsigned& j)
2015 {
2016 return Haloed_node_pt[p][j];
2017 }
2018
2019 /// Add haloed node whose halo counterpart is held
2020 /// on processor p to the storage scheme for haloed nodes.
2021 void add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
2022 {
2023 Haloed_node_pt[p].push_back(nod_pt);
2024 }
2025
2026 /// Bool for output of halo elements
2028
2029
2030 /// Function to suppress resizing of halo nodes -- optmisation
2031 /// but call it at your own risk!
2033 {
2035 }
2036
2037
2038 /// Function to (re-)enable resizing of halo nodes -- this returns
2039 /// things to the default behaviour.
2041 {
2043 }
2044
2045 /// Function to disable halo element output
2047 {
2048 Output_halo_elements = false;
2049 }
2050
2051 /// Function to enable halo element output
2053 {
2054 Output_halo_elements = true;
2055 }
2056
2057 /// Total number of shared nodes in this Mesh
2058 unsigned nshared_node()
2059 {
2060 unsigned n = 0;
2061 for (std::map<unsigned, Vector<Node*>>::iterator it =
2062 Shared_node_pt.begin();
2063 it != Shared_node_pt.end();
2064 it++)
2065 {
2066 n += it->second.size();
2067 }
2068 return n;
2069 }
2070
2071 /// Doc shared nodes
2073 {
2074 for (std::map<unsigned, Vector<Node*>>::iterator it =
2075 Shared_node_pt.begin();
2076 it != Shared_node_pt.end();
2077 it++)
2078 {
2079 unsigned n = it->second.size();
2080 for (unsigned j = 0; j < n; j++)
2081 {
2082 oomph_info << "Shared node with proc " << it->first << " ";
2083 Node* nod_pt = it->second[j];
2084 unsigned ndim = nod_pt->ndim();
2085 for (unsigned i = 0; i < ndim; i++)
2086 {
2087 oomph_info << nod_pt->x(i) << " ";
2088 }
2089 oomph_info << nod_pt->is_hanging() << " ";
2090 oomph_info << j << " " << nod_pt << std::endl;
2091 }
2092 }
2093 }
2094
2095 /// Number of shared nodes in this Mesh who have a counterpart
2096 /// on processor p.
2097 unsigned nshared_node(const unsigned& p)
2098 {
2099 // Memory saving version of: return Shared_node_pt[p].size();
2100 std::map<unsigned, Vector<Node*>>::iterator it = Shared_node_pt.find(p);
2101 if (it == Shared_node_pt.end())
2102 {
2103 return 0;
2104 }
2105 return (*it).second.size();
2106 }
2107
2108 /// Access fct to the j-th shared node in this Mesh
2109 /// who has a counterpart on processor p.
2110 Node* shared_node_pt(const unsigned& p, const unsigned& j)
2111 {
2112 return Shared_node_pt[p][j];
2113 }
2114
2115 /// Get vector of pointers to shared nodes with processor p.
2116 /// Required for faster search in
2117 /// Missing_masters_functions::add_external_haloed_node_helper() and
2118 /// Missing_masters_functions::add_external_haloed_master_node_helper()
2120 {
2121 unsigned np = nshared_node(p);
2123 for (unsigned j = 0; j < np; j++)
2124 {
2125 shared_node_pt[j] = Shared_node_pt[p][j];
2126 }
2127 }
2128
2129
2130 /// Add shared node whose counterpart is held
2131 /// on processor p to the storage scheme for shared nodes.
2132 /// (NB: ensure that this routine is called twice, once for each process)
2133 void add_shared_node_pt(const unsigned& p, Node*& nod_pt)
2134 {
2135 Shared_node_pt[p].push_back(nod_pt);
2136 }
2137
2138 /// Get halo node stats for this distributed mesh:
2139 /// Average/max/min number of halo nodes over all processors.
2140 /// \b Careful: Involves MPI Broadcasts and must therefore
2141 /// be called on all processors!
2142 void get_halo_node_stats(double& av_number,
2143 unsigned& max_number,
2144 unsigned& min_number);
2145
2146 /// Get haloed node stats for this distributed mesh:
2147 /// Average/max/min number of haloed nodes over all processors.
2148 /// \b Careful: Involves MPI Broadcasts and must therefore
2149 /// be called on all processors!
2150 void get_haloed_node_stats(double& av_number,
2151 unsigned& max_number,
2152 unsigned& min_number);
2153
2154 // External halo(ed) elements are "source/other" elements which are
2155 // on different processes to the element for which they are the source
2156
2157 /// Output all external halo elements
2158 void output_external_halo_elements(std::ostream& outfile,
2159 const unsigned& n_plot = 5)
2160 {
2161 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2163 it != External_halo_element_pt.end();
2164 it++)
2165 {
2166 unsigned p = (*it).first;
2167 output_external_halo_elements(p, outfile, n_plot);
2168 }
2169 }
2170
2171 /// Output all external halo elements with processor p
2172 void output_external_halo_elements(const unsigned& p,
2173 std::ostream& outfile,
2174 const unsigned& n_plot = 5)
2175 {
2176 unsigned nel = External_halo_element_pt[p].size();
2177 for (unsigned e = 0; e < nel; e++)
2178 {
2179 FiniteElement* fe_pt =
2180 dynamic_cast<FiniteElement*>(External_halo_element_pt[p][e]);
2181 if (fe_pt != 0)
2182 {
2183 fe_pt->output(outfile, n_plot);
2184 }
2185 }
2186 }
2187
2188
2189 /// Output all external haloed elements
2190 void output_external_haloed_elements(std::ostream& outfile,
2191 const unsigned& n_plot = 5)
2192 {
2193 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2195 it != External_haloed_element_pt.end();
2196 it++)
2197 {
2198 unsigned p = (*it).first;
2199 output_external_haloed_elements(p, outfile, n_plot);
2200 }
2201 }
2202
2203 /// Output all external haloed elements with processor p
2204 void output_external_haloed_elements(const unsigned& p,
2205 std::ostream& outfile,
2206 const unsigned& n_plot = 5)
2207 {
2208 unsigned nel = External_haloed_element_pt[p].size();
2209 for (unsigned e = 0; e < nel; e++)
2210 {
2211 FiniteElement* fe_pt =
2212 dynamic_cast<FiniteElement*>(External_haloed_element_pt[p][e]);
2213 if (fe_pt != 0)
2214 {
2215 fe_pt->output(outfile, n_plot);
2216 }
2217 }
2218 }
2219
2220
2221 /// Total number of external halo elements in this Mesh
2223 {
2224 unsigned n = 0;
2225 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2227 it != External_halo_element_pt.end();
2228 it++)
2229 {
2230 n += it->second.size();
2231 }
2232 return n;
2233 }
2234
2235 /// Number of external halo elements in this Mesh whose non-halo
2236 /// counterpart is held on processor p.
2237 unsigned nexternal_halo_element(const unsigned& p)
2238 {
2239 // Memory saving version of: return External_halo_element_pt[p].size();
2240 std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2242 if (it == External_halo_element_pt.end())
2243 {
2244 return 0;
2245 }
2246 return (*it).second.size();
2247 }
2248
2249 /// Access fct to the e-th external halo element in this Mesh
2250 /// whose non-halo counterpart is held on processor p.
2252 const unsigned& e)
2253 {
2254 return External_halo_element_pt[p][e];
2255 }
2256
2257 /// Add external halo element whose non-halo counterpart is held
2258 /// on processor p to this Mesh.
2259 void add_external_halo_element_pt(const unsigned& p,
2260 GeneralisedElement*& el_pt)
2261 {
2262 External_halo_element_pt[p].push_back(el_pt);
2263 el_pt->set_halo(p);
2264 }
2265
2266 /// Total number of external haloed elements in this Mesh
2268 {
2269 unsigned n = 0;
2270 for (std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2272 it != External_haloed_element_pt.end();
2273 it++)
2274 {
2275 n += it->second.size();
2276 }
2277 return n;
2278 }
2279
2280 /// Number of external haloed elements in this Mesh whose non-halo
2281 /// counterpart is held on processor p.
2282 unsigned nexternal_haloed_element(const unsigned& p)
2283 {
2284 // Memory saving version of: return External_haloed_element_pt[p].size();
2285 std::map<unsigned, Vector<GeneralisedElement*>>::iterator it =
2287 if (it == External_haloed_element_pt.end())
2288 {
2289 return 0;
2290 }
2291 return (*it).second.size();
2292 }
2293
2294 /// Access fct to the e-th external haloed element in this Mesh
2295 /// whose non-halo counterpart is held on processor p.
2297 const unsigned& e)
2298 {
2299 return External_haloed_element_pt[p][e];
2300 }
2301
2302 /// Add external haloed element whose non-halo counterpart is held
2303 /// on processor p to the storage scheme for haloed elements.
2304 unsigned add_external_haloed_element_pt(const unsigned& p,
2305 GeneralisedElement*& el_pt);
2306
2307 /// Total number of external halo nodes in this Mesh
2309 {
2310 unsigned n = 0;
2311 for (std::map<unsigned, Vector<Node*>>::iterator it =
2312 External_halo_node_pt.begin();
2313 it != External_halo_node_pt.end();
2314 it++)
2315 {
2316 n += it->second.size();
2317 }
2318 return n;
2319 }
2320
2321
2322 /// Get vector of pointers to all external halo nodes
2324 {
2325 unsigned n_total = nexternal_halo_node();
2326 external_halo_node_pt.reserve(n_total);
2327 for (std::map<unsigned, Vector<Node*>>::iterator it =
2328 External_halo_node_pt.begin();
2329 it != External_halo_node_pt.end();
2330 it++)
2331 {
2332 unsigned np = (it->second).size();
2333 for (unsigned j = 0; j < np; j++)
2334 {
2335 external_halo_node_pt.push_back((it->second)[j]);
2336 }
2337 }
2338#ifdef PARANOID
2339 if (external_halo_node_pt.size() != n_total)
2340 {
2341 std::ostringstream error_stream;
2342 error_stream << "Total number of external halo nodes, " << n_total
2343 << " doesn't match number of entries \n in vector, "
2344 << external_halo_node_pt.size() << std::endl;
2345 throw OomphLibError(
2346 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2347 }
2348#endif
2349 }
2350
2351
2352 /// Number of external halo nodes in this Mesh whose non-halo
2353 /// (external) counterpart is held on processor p.
2354 unsigned nexternal_halo_node(const unsigned& p)
2355 {
2356 // Memory saving version of: return External_halo_node_pt[p].size();
2357 std::map<unsigned, Vector<Node*>>::iterator it =
2358 External_halo_node_pt.find(p);
2359 if (it == External_halo_node_pt.end())
2360 {
2361 return 0;
2362 }
2363 return (*it).second.size();
2364 }
2365
2366 /// Add external halo node whose non-halo (external) counterpart
2367 /// is held on processor p to the storage scheme for halo nodes.
2368 void add_external_halo_node_pt(const unsigned& p, Node*& nod_pt)
2369 {
2370 External_halo_node_pt[p].push_back(nod_pt);
2371 nod_pt->set_halo(p);
2372 }
2373
2374
2375 /// Access fct to the j-th external halo node in this Mesh
2376 /// whose non-halo external counterpart is held on processor p.
2377 Node*& external_halo_node_pt(const unsigned& p, const unsigned& j)
2378 {
2379 return External_halo_node_pt[p][j];
2380 }
2381
2382 /// Access fct to vector of external halo node in this Mesh
2383 /// whose non-halo external counterpart is held on processor p. (read only)
2385 {
2386 std::map<unsigned, Vector<Node*>>::iterator it =
2387 External_halo_node_pt.find(p);
2388 if (it == External_halo_node_pt.end())
2389 {
2390 Vector<Node*> tmp;
2391 return tmp;
2392 }
2393 return (*it).second;
2394 }
2395
2396 /// Set vector of external halo node in this Mesh
2397 /// whose non-halo external counterpart is held on processor p.
2398 void set_external_halo_node_pt(const unsigned& p,
2400 {
2402 }
2403
2404 /// Null out specified external halo node (used when deleting duplicates)
2405 void null_external_halo_node(const unsigned& p, Node* nod_pt);
2406
2407 /// Consolidate external halo node storage by removing nulled out
2408 /// pointes in external halo and haloed schemes
2410
2411 /// Total number of external haloed nodes in this Mesh
2413 {
2414 unsigned n = 0;
2415 for (std::map<unsigned, Vector<Node*>>::iterator it =
2417 it != External_haloed_node_pt.end();
2418 it++)
2419 {
2420 n += it->second.size();
2421 }
2422 return n;
2423 }
2424
2425 /// Number of external haloed nodes in this Mesh
2426 /// whose halo (external) counterpart is held on processor p.
2427 unsigned nexternal_haloed_node(const unsigned& p)
2428 {
2429 // Memory saving version of: return External_haloed_node_pt[p].size();
2430 std::map<unsigned, Vector<Node*>>::iterator it =
2432 if (it == External_haloed_node_pt.end())
2433 {
2434 return 0;
2435 }
2436 return (*it).second.size();
2437 }
2438
2439 /// Access fct to the j-th external haloed node in this Mesh
2440 /// whose halo external counterpart is held on processor p.
2441 Node*& external_haloed_node_pt(const unsigned& p, const unsigned& j)
2442 {
2443 return External_haloed_node_pt[p][j];
2444 }
2445
2446 /// Add external haloed node whose halo (external) counterpart
2447 /// is held on processor p to the storage scheme for haloed nodes.
2448 unsigned add_external_haloed_node_pt(const unsigned& p, Node*& nod_pt);
2449
2450 /// Access fct to vector of external haloed node in this Mesh
2451 /// whose halo external counterpart is held on processor p. (read only)
2453 {
2454 std::map<unsigned, Vector<Node*>>::iterator it =
2456 if (it == External_haloed_node_pt.end())
2457 {
2458 Vector<Node*> tmp;
2459 return tmp;
2460 }
2461 return (*it).second;
2462 }
2463
2464 /// Set vector of external haloed node in this Mesh
2465 /// whose halo external counterpart is held on processor p.
2467 const unsigned& p, const Vector<Node*>& external_haloed_node_pt)
2468 {
2470 }
2471
2472 /// Return the set of processors that hold external halo nodes. This
2473 /// is required to avoid having to pass a communicator into the node_update
2474 /// functions for Algebraic-based and MacroElement-based Meshes
2475 std::set<int> external_halo_proc()
2476 {
2477 std::set<int> procs;
2478 for (std::map<unsigned, Vector<Node*>>::iterator it =
2479 External_halo_node_pt.begin();
2480 it != External_halo_node_pt.end();
2481 it++)
2482 {
2483 procs.insert((*it).first);
2484 }
2485 return procs;
2486 }
2487
2488 /// Creates the shared boundaries, only used in unstructured meshes
2489 /// In this case with the "TriangleMesh" class
2491 OomphCommunicator* comm_pt,
2492 const Vector<unsigned>& element_domain,
2493 const Vector<GeneralisedElement*>& backed_up_el_pt,
2494 const Vector<FiniteElement*>& backed_up_f_el_pt,
2495 std::map<Data*, std::set<unsigned>>& processors_associated_with_data,
2496 const bool& overrule_keep_as_halo_element_status)
2497 {
2498 std::ostringstream error_stream;
2499 error_stream << "Empty default create_shared_boundaries() method"
2500 << "called.\n";
2501 error_stream << "This should be overloaded in a specific "
2502 << "TriangleMeshBase\n";
2503 throw OomphLibError(error_stream.str(),
2504 "Mesh::create_shared_boundaries()",
2505 OOMPH_EXCEPTION_LOCATION);
2506 }
2507
2508 // Check if necessary to add the element as haloed or if it has been
2509 // previously added to the haloed scheme
2511 const unsigned& p, GeneralisedElement*& el_pt)
2512 {
2513 std::ostringstream error_stream;
2514 error_stream << "Empty default try_to_add_root_haloed_element_pt() method"
2515 << "called.\n";
2516 error_stream << "This should be overloaded in a specific "
2517 << "TriangleMeshBase\n";
2518 throw OomphLibError(error_stream.str(),
2519 "Mesh::try_to_add_root_haloed_element_pt()",
2520 OOMPH_EXCEPTION_LOCATION);
2521 }
2522
2523 // Check if necessary to add the node as haloed or if it has been
2524 // previously added to the haloed scheme
2525 virtual unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
2526 {
2527 std::ostringstream error_stream;
2528 error_stream << "Empty default try_to_add_haloed_node_pt() method"
2529 << "called.\n";
2530 error_stream << "This should be overloaded in a specific "
2531 << "TriangleMeshBase\n";
2532 throw OomphLibError(error_stream.str(),
2533 "Mesh::try_to_add_haloed_node_pt()",
2534 OOMPH_EXCEPTION_LOCATION);
2535 }
2536
2537#endif
2538
2539 /// Wipe the storage for all externally-based elements
2541 };
2542
2543
2544 /// ///////////////////////////////////////////////////////////////////////
2545 /// ///////////////////////////////////////////////////////////////////////
2546 /// ///////////////////////////////////////////////////////////////////////
2547
2548 class SolidICProblem;
2549
2550 //========================================================================
2551 /// General SolidMesh class.
2552 ///
2553 /// Solid meshes contain SolidFiniteElements which contain
2554 /// SolidNodes. This class simply overloads the appropriate access
2555 /// functions from the underlying Mesh class.
2556 //
2557 // Needs to be derived from Mesh with virtual so that
2558 // solid meshes can be derived from general meshes, without
2559 // multiple copies of Mesh objects
2560 //========================================================================
2561 class SolidMesh : public virtual Mesh
2562 {
2563 public:
2564 /// Default constructor
2566
2567 /// Broken copy constructor
2568 SolidMesh(const SolidMesh& dummy) = delete;
2569
2570 /// Broken assignment operator
2571 void operator=(const SolidMesh&) = delete;
2572
2573 /// Constructor builds combined mesh from the meshes specified.
2574 /// Note: This simply merges the meshes' elements and nodes (ignoring
2575 /// duplicates; no boundary information etc. is created).
2576 SolidMesh(const Vector<SolidMesh*>& sub_mesh_pt)
2577 {
2578#ifdef OOMPH_HAS_MPI
2579 // Mesh hasn't been distributed: Null out pointer to communicator
2580 Comm_pt = 0;
2581#endif
2582
2583 unsigned n = sub_mesh_pt.size();
2584 Vector<Mesh*> sub_mesh_mesh_pt(n);
2585 for (unsigned i = 0; i < n; i++)
2586 {
2587 sub_mesh_mesh_pt[i] = static_cast<Mesh*>(sub_mesh_pt[i]);
2588 }
2589 merge_meshes(sub_mesh_mesh_pt);
2590 }
2591
2592 /// Return a pointer to the n-th global SolidNode
2593 // Can safely cast the nodes to SolidNodes
2594 SolidNode* node_pt(const unsigned long& n)
2595 {
2596#ifdef PARANOID
2597 if (!dynamic_cast<SolidNode*>(Node_pt[n]))
2598 {
2599 std::ostringstream error_stream;
2600 error_stream << "Error: Node " << n << "is a "
2601 << typeid(Node_pt[n]).name() << ", not an SolidNode"
2602 << std::endl;
2603 throw OomphLibError(
2604 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2605 }
2606#endif
2607 // Return a static cast to the Node_pt
2608 return (static_cast<SolidNode*>(Node_pt[n]));
2609 }
2610
2611 /// Return n-th SolidNodes on b-th boundary
2612 SolidNode* boundary_node_pt(const unsigned& b, const unsigned& n)
2613 {
2614#ifdef PARANOID
2615 if (!dynamic_cast<SolidNode*>(Mesh::boundary_node_pt(b, n)))
2616 {
2617 std::ostringstream error_stream;
2618 error_stream << "Error: Node " << n << "is a "
2619 << typeid(Mesh::boundary_node_pt(b, n)).name()
2620 << ", not an SolidNode" << std::endl;
2621 throw OomphLibError(
2622 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2623 }
2624#endif
2625 return static_cast<SolidNode*>(Mesh::boundary_node_pt(b, n));
2626 }
2627
2628 /// Return the n-th local SolidNode in elemnet e.
2629 /// This is required to cast the nodes in a solid mesh to be
2630 /// SolidNodes and therefore allow access to the extra SolidNode data
2631 SolidNode* element_node_pt(const unsigned long& e, const unsigned& n)
2632 {
2633#ifdef PARANOID
2634 // Try to cast to FiniteElement
2635 FiniteElement* el_pt = dynamic_cast<FiniteElement*>(Element_pt[e]);
2636 if (el_pt == 0)
2637 {
2638 // Error
2639 throw OomphLibError("Failed cast to FiniteElement* ",
2640 OOMPH_CURRENT_FUNCTION,
2641 OOMPH_EXCEPTION_LOCATION);
2642 }
2643 if (!dynamic_cast<SolidNode*>(el_pt->node_pt(n)))
2644 {
2645 std::ostringstream error_message;
2646 Node* np = el_pt->node_pt(n);
2647 error_message << "Error: Node " << n << " of element " << e << "is a "
2648 << typeid(*np).name() << ", not an SolidNode"
2649 << std::endl;
2650
2651 throw OomphLibError(error_message.str(),
2652 OOMPH_CURRENT_FUNCTION,
2653 OOMPH_EXCEPTION_LOCATION);
2654 }
2655#endif
2656 // Return a cast to an SolidNode
2657 return (static_cast<SolidNode*>(
2658 dynamic_cast<FiniteElement*>(Element_pt[e])->node_pt(n)));
2659 }
2660
2661 /// Make the current configuration the undeformed one by
2662 /// setting the nodal Lagrangian coordinates to their current
2663 /// Eulerian ones
2665
2666
2667 /// Scale all nodal coordinates by given factor and re-assign the
2668 /// Lagrangian coordinates
2669 void scale_mesh(const double& factor)
2670 {
2671 Mesh::scale_mesh(factor);
2672
2673 // Re-assign the Lagrangian coordinates
2675 }
2676 /// Static problem that can be used to assign initial conditions
2677 /// on a given solid mesh (need to define this as a static problem
2678 /// somewhere because deleting the problem would wipe out the mesh too!)
2680 };
2681
2682
2683 /// ////////////////////////////////////////////////////////////////////////
2684 /// ////////////////////////////////////////////////////////////////////////
2685 /// ////////////////////////////////////////////////////////////////////////
2686
2687
2688 //===================================================
2689 /// Edge class
2690 //===================================================
2691 class Edge
2692 {
2693 public:
2694 /// Constructor: Pass in the two vertex nodes
2696 {
2697 if (node1_pt == node2_pt)
2698 {
2699#ifdef PARANOID
2700 std::ostringstream error_stream;
2701 error_stream << "Edge cannot have two identical vertex nodes\n";
2702 throw OomphLibError(
2703 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2704#endif
2705 }
2706
2707
2708 // Sort lexicographically based on pointer address of nodes
2709 if (node1_pt > node2_pt)
2710 {
2713 }
2714 else
2715 {
2718 }
2719 }
2720
2721
2722 /// Access to the first vertex node
2724 {
2725 return Node1_pt;
2726 }
2727
2728 /// Access to the second vertex node
2730 {
2731 return Node2_pt;
2732 }
2733
2734 /// Comparison operator
2735 bool operator==(const Edge& other) const
2736 {
2737 if ((Node1_pt == other.node1_pt()) && (Node2_pt == other.node2_pt()))
2738
2739 {
2740 return true;
2741 }
2742 else
2743 {
2744 return false;
2745 }
2746 }
2747
2748
2749 /// Less-than operator
2750 bool operator<(const Edge& other) const
2751 {
2752 if (Node1_pt < other.node1_pt())
2753 {
2754 return true;
2755 }
2756 else if (Node1_pt == other.node1_pt())
2757 {
2758 if (Node2_pt < other.node2_pt())
2759 {
2760 return true;
2761 }
2762 else
2763 {
2764 return false;
2765 }
2766 }
2767 else
2768 {
2769 return false;
2770 }
2771 }
2772
2773 /// Test whether the Edge lies on a boundary. Relatively simple
2774 /// test, based on both vertices lying on (some) boundary.
2775 bool is_on_boundary() const
2776 {
2778 }
2779
2780
2781 /// Test whether the Edge is a boundary edge, i.e. does it
2782 /// connnect two boundary nodes?
2783 bool is_boundary_edge() const
2784 {
2785 return ((dynamic_cast<BoundaryNodeBase*>(Node1_pt) != 0) &&
2786 (dynamic_cast<BoundaryNodeBase*>(Node2_pt) != 0));
2787 }
2788
2789 private:
2790 /// First vertex node
2792
2793 /// Second vertex node
2795 };
2796
2797
2798 /// ////////////////////////////////////////////////////////////////////
2799 /// ////////////////////////////////////////////////////////////////////
2800 /// ////////////////////////////////////////////////////////////////////
2801
2802
2803 //=================================================================
2804 /// Namespace with helper function to check element type in mesh
2805 /// constructors (say).
2806 //=================================================================
2807 namespace MeshChecker
2808 {
2809 //=================================================================
2810 /// Helper function to assert that finite element of type ELEMENT
2811 /// can be cast to base class of type GEOM_ELEMENT_BASE -- useful
2812 /// to avoid confusion if a mesh that was written for a specific
2813 /// element type (e.g. a QElement) is used with another one (e.g.
2814 /// a TElement. First argument specifies the required spatial dimension
2815 /// of the element (i.e. the number of local coordinates). The optional
2816 /// second argument specifies the required nnode_1d (i.e. the number
2817 /// of nodes along a 1D element edge). Can be omitted if the mesh
2818 /// can handle any number in which case this test is skipped.
2819 //=================================================================
2820 template<class GEOM_ELEMENT_BASE, class ELEMENT>
2821 void assert_geometric_element(const unsigned& dim,
2822 const unsigned& nnode_1d = 0)
2823 {
2824 // Only do tests in paranoia mode
2825#ifndef PARANOID
2826 return;
2827#endif
2828
2829 // Instantiate element
2830 ELEMENT* el_pt = new ELEMENT;
2831
2832 // Can we cast to required geometric element base type
2833 if (dynamic_cast<GEOM_ELEMENT_BASE*>(el_pt) == 0)
2834 {
2835 std::stringstream error_message;
2836 error_message << "You have specified an illegal element type! Element "
2837 "is of type \n\n"
2838 << typeid(el_pt).name()
2839 << "\n\nand cannot be cast to type \n\n "
2840 << typeid(GEOM_ELEMENT_BASE).name() << "\n\n";
2841 throw OomphLibError(error_message.str(),
2842 OOMPH_CURRENT_FUNCTION,
2843 OOMPH_EXCEPTION_LOCATION);
2844 }
2845
2846 // Does the dimension match?
2847 if (dim != el_pt->dim())
2848 {
2849 std::stringstream error_message;
2850 error_message << "You have specified an illegal element type! Element "
2851 "is of type \n\n"
2852 << typeid(el_pt).name()
2853 << "\n\nand has dimension = " << el_pt->dim()
2854 << " but we need dim = " << dim << std::endl;
2855 throw OomphLibError(error_message.str(),
2856 OOMPH_CURRENT_FUNCTION,
2857 OOMPH_EXCEPTION_LOCATION);
2858 }
2859
2860 // Does nnode_1d match?
2861 if (nnode_1d != 0)
2862 {
2863 if (nnode_1d != el_pt->nnode_1d())
2864 {
2865 std::stringstream error_message;
2866 error_message << "You have specified an illegal element type! "
2867 "Element is of type \n\n"
2868 << typeid(el_pt).name()
2869 << "\n\nand has nnode_1d = " << el_pt->nnode_1d()
2870 << " but we need nnode_1d = " << nnode_1d << std::endl;
2871 throw OomphLibError(error_message.str(),
2872 OOMPH_CURRENT_FUNCTION,
2873 OOMPH_EXCEPTION_LOCATION);
2874 }
2875 }
2876
2877 // Clean up
2878 delete el_pt;
2879 }
2880
2881 } // namespace MeshChecker
2882
2883
2884 /// ////////////////////////////////////////////////////////////////////
2885 /// ////////////////////////////////////////////////////////////////////
2886 /// ////////////////////////////////////////////////////////////////////
2887
2888
2889 //=================paraview_helper=================================
2890 /// Namespace for paraview-style output helper functions
2891 //=================================================================
2892 namespace ParaviewHelper
2893 {
2894 /// Write the pvd file header
2895 extern void write_pvd_header(std::ofstream& pvd_file);
2896
2897 /// Add name of output file and associated continuous time
2898 /// to pvd file.
2899 extern void write_pvd_information(std::ofstream& pvd_file,
2900 const std::string& output_filename,
2901 const double& time);
2902
2903 /// Write the pvd file footer
2904 extern void write_pvd_footer(std::ofstream& pvd_file);
2905
2906 } // namespace ParaviewHelper
2907
2908 /// /////////////////////////////////////////////////////////////
2909 /// /////////////////////////////////////////////////////////////
2910 /// /////////////////////////////////////////////////////////////
2911
2912 namespace NodeOrdering
2913 {
2914 /// Function for ordering nodes. Return true if first node's position is
2915 /// "before" second nodes. Dimension 0 checked first, then... until they
2916 /// are different (by more than tol=1e-10). If they are both in exactly
2917 /// the same place an error is thrown.
2918 inline bool node_global_position_comparison(Node* nd1_pt, Node* nd2_pt)
2919 {
2920 unsigned ndim = nd1_pt->ndim();
2921
2922 unsigned j;
2923 for (j = 0; j < ndim; j++)
2924 {
2925 if (std::abs(nd1_pt->x(j) - nd2_pt->x(j)) > 1e-10)
2926 {
2927 if (nd1_pt->x(j) < nd2_pt->x(j))
2928 {
2929 return true;
2930 }
2931 else
2932 {
2933 return false;
2934 }
2935 }
2936 // otherwise they are the same in this dimension, check next one.
2937 }
2938
2939 std::string err = "Nodes are at the same point to ~ 1e-10!";
2940 err += " difference is " +
2941 StringConversion::to_string(std::abs(nd1_pt->x(j) - nd2_pt->x(j)));
2942 throw OomphLibError(
2943 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2944 }
2945 } // namespace NodeOrdering
2946
2947
2948} // namespace oomph
2949
2950#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
GeneralisedTimestepper used to store the arclength derivatives and pervious solutions required in con...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
virtual bool is_a_copy() const
Return a boolean to indicate whether the Data objact contains any copied values. A base Data object c...
Definition: nodes.h:253
void set_halo(const unsigned &non_halo_proc_ID)
Label the node as halo and specify processor that holds non-halo counterpart.
Definition: nodes.h:520
Information for documentation of results: Directory and file number to enable output in the form RESL...
void disable_doc()
Disable documentation.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5014
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
Definition: mesh.h:2692
Node * node2_pt() const
Access to the second vertex node.
Definition: mesh.h:2729
bool is_on_boundary() const
Test whether the Edge lies on a boundary. Relatively simple test, based on both vertices lying on (so...
Definition: mesh.h:2775
Node * node1_pt() const
Access to the first vertex node.
Definition: mesh.h:2723
bool is_boundary_edge() const
Test whether the Edge is a boundary edge, i.e. does it connnect two boundary nodes?
Definition: mesh.h:2783
Node * Node1_pt
First vertex node.
Definition: mesh.h:2791
Edge(Node *node1_pt, Node *node2_pt)
Constructor: Pass in the two vertex nodes.
Definition: mesh.h:2695
bool operator<(const Edge &other) const
Less-than operator.
Definition: mesh.h:2750
Node * Node2_pt
Second vertex node.
Definition: mesh.h:2794
bool operator==(const Edge &other) const
Comparison operator.
Definition: mesh.h:2735
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4482
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
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4675
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
virtual void compute_error(FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Calculate the norm of the error and that of the exact solution.
Definition: elements.h:3198
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
A Generalised Element class.
Definition: elements.h:73
bool is_halo() const
Is this element a halo?
Definition: elements.h:1163
void set_halo(const unsigned &non_halo_proc_ID)
Label the element as halo and specify processor that holds non-halo counterpart.
Definition: elements.h:1151
virtual void compute_norm(Vector< double > &norm)
Compute norm of solution – broken virtual can be overloaded by element writer to implement whatever n...
Definition: elements.h:1128
A general mesh class.
Definition: mesh.h:67
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
void set_external_halo_node_pt(const unsigned &p, const Vector< Node * > &external_halo_node_pt)
Set vector of external halo node in this Mesh whose non-halo external counterpart is held on processo...
Definition: mesh.h:2398
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
void(FiniteElement::* SteadyExactSolutionFctPt)(const Vector< double > &x, Vector< double > &soln)
Typedef for function pointer to function that computes steady exact solution.
Definition: mesh.h:223
Vector< Node * > external_haloed_node_pt(const unsigned &p)
Access fct to vector of external haloed node in this Mesh whose halo external counterpart is held on ...
Definition: mesh.h:2452
GeneralisedElement * element_pt(const unsigned long &e) const
Return pointer to element e (const version)
Definition: mesh.h:454
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
void enable_resizing_of_halo_nodes()
Function to (re-)enable resizing of halo nodes – this returns things to the default behaviour.
Definition: mesh.h:2040
std::map< unsigned, Vector< Node * > > Shared_node_pt
Map of vectors holding the pointers to the shared nodes. These are all the nodes that are on two "nei...
Definition: mesh.h:116
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
bool does_pointer_correspond_to_mesh_data(double *const &parameter_pt)
Does the double pointer correspond to any mesh data.
Definition: mesh.cc:2471
virtual void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1351
void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
Distribute the problem Add to vector of pointers to deleted elements.
Definition: mesh.h:1645
void disable_resizing_of_halo_nodes()
Function to suppress resizing of halo nodes – optmisation but call it at your own risk!
Definition: mesh.h:2032
unsigned ndof_types() const
Return number of dof types in mesh.
Definition: mesh.cc:8764
void remove_boundary_node(const unsigned &b, Node *const &node_pt)
Remove a node from the boundary b.
Definition: mesh.cc:221
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
void output_external_halo_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external halo elements.
Definition: mesh.h:2158
void setup_shared_node_scheme()
Setup shared node scheme.
Definition: mesh.cc:2712
Node *& external_halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external halo node in this Mesh whose non-halo external counterpart is held on...
Definition: mesh.h:2377
GeneralisedElement *& external_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external halo element in this Mesh whose non-halo counterpart is held on proce...
Definition: mesh.h:2251
void max_and_min_element_size(double &max_size, double &min_size)
Determine max and min area for all FiniteElements in the mesh (non-FiniteElements are ignored)
Definition: mesh.h:692
std::map< unsigned, Vector< GeneralisedElement * > > Root_haloed_element_pt
Map of vectors holding the pointers to the root haloed elements.
Definition: mesh.h:103
void set_keep_all_elements_as_halos()
Call this function to keep all the elements as halo elements.
Definition: mesh.h:1621
unsigned nroot_haloed_element(const unsigned &p)
Number of root haloed elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1936
virtual unsigned try_to_add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Definition: mesh.h:2525
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition: mesh.h:91
void add_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add halo node whose non-halo counterpart is held on processor p to the storage scheme for halo nodes.
Definition: mesh.h:1906
unsigned nexternal_haloed_element()
Total number of external haloed elements in this Mesh.
Definition: mesh.h:2267
unsigned add_external_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add external haloed node whose halo (external) counterpart is held on processor p to the storage sche...
Definition: mesh.cc:9516
unsigned nexternal_halo_node()
Total number of external halo nodes in this Mesh.
Definition: mesh.h:2308
GeneralisedElement *& root_halo_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th root halo element in this Mesh whose non-halo counterpart is held on processor...
Definition: mesh.h:1862
unsigned nexternal_haloed_element(const unsigned &p)
Number of external haloed elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:2282
unsigned nodal_dimension() const
Return number of nodal dimension in mesh.
Definition: mesh.cc:9055
void add_external_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external halo element whose non-halo counterpart is held on processor p to this Mesh.
Definition: mesh.h:2259
virtual void compute_norm(double &norm)
Compute norm of solution by summing contributions of compute_norm(...) for all constituent elements i...
Definition: mesh.h:1068
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned > > &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Virtual function to perform the reset boundary elements info rutines.
Definition: mesh.h:284
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
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them....
Definition: mesh.h:407
Node * boundary_node_pt(const unsigned &b, const unsigned &n) const
Return pointer to node n on boundary b.
Definition: mesh.h:499
std::map< unsigned, Vector< GeneralisedElement * > > External_haloed_element_pt
Map of vectors holding the pointers to the external haloed elements.
Definition: mesh.h:129
void add_shared_node_pt(const unsigned &p, Node *&nod_pt)
Add shared node whose counterpart is held on processor p to the storage scheme for shared nodes....
Definition: mesh.h:2133
void flush_element_storage()
Flush storage for elements (only) by emptying the vectors that store the pointers to them....
Definition: mesh.h:423
void output(const std::string &output_filename)
Output for all elements.
Definition: mesh.h:975
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
virtual void setup_boundary_element_info(std::ostream &outfile)
Setup lookup schemes which establish whic elements are located next to mesh's boundaries....
Definition: mesh.h:281
GeneralisedElement *& external_haloed_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th external haloed element in this Mesh whose non-halo counterpart is held on pro...
Definition: mesh.h:2296
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
void get_haloed_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get haloed node stats for this distributed mesh: Average/max/min number of haloed nodes over all proc...
Definition: mesh.cc:4904
void(FiniteElement::* UnsteadyExactSolutionFctPt)(const double &time, const Vector< double > &x, Vector< double > &soln)
Typedef for function pointer to function that computes unsteady exact solution.
Definition: mesh.h:228
void check_halo_schemes(DocInfo &doc_info, double &max_permitted_error_for_halo_check)
Check halo and shared schemes on the mesh.
Definition: mesh.cc:6881
virtual void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Function that can be used to set any additional timestepper data stored at the Mesh (as opposed to no...
Definition: mesh.cc:2402
unsigned nhaloed_node(const unsigned &p)
Number of haloed nodes in this Mesh whose haloed counterpart is held on processor p.
Definition: mesh.h:2001
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the elements. The ostream specifies the output stream to which...
Definition: mesh.cc:746
virtual void classify_halo_and_haloed_nodes(DocInfo &doc_info, const bool &report_stats)
Classify the halo and haloed nodes in the mesh. Virtual so it can be overloaded to perform additional...
Definition: mesh.cc:3262
void add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root haloed element whose non-halo counterpart is held on processor p to the storage scheme for h...
Definition: mesh.h:1977
void synchronise_shared_nodes(const bool &report_stats)
Synchronise shared node lookup schemes to cater for the the case where: (1) a certain node on the cur...
Definition: mesh.cc:2890
void check_inverted_elements(bool &mesh_has_inverted_elements)
Check for inverted elements and report outcome in boolean variable. This visits all elements at their...
Definition: mesh.h:742
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
void output_external_halo_elements(const unsigned &p, std::ostream &outfile, const unsigned &n_plot=5)
Output all external halo elements with processor p.
Definition: mesh.h:2172
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
unsigned nexternal_haloed_node(const unsigned &p)
Number of external haloed nodes in this Mesh whose halo (external) counterpart is held on processor p...
Definition: mesh.h:2427
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
void copy_boundary_node_data_from_nodes()
Replace existing boundary node lookup schemes with new schemes created using the boundary data stored...
Definition: mesh.h:526
void set_elemental_internal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the internal data stored within elements in the meah.
Definition: mesh.cc:2535
void doc_boundary_coordinates(const unsigned &b, std::ofstream &the_file)
Output boundary coordinates on boundary b – template argument specifies the bulk element type (needed...
Definition: mesh.h:303
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
bool Keep_all_elements_as_halos
bool to indicate whether to keep all elements in a mesh as halos or not
Definition: mesh.h:141
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
std::map< unsigned, Vector< Node * > > External_haloed_node_pt
Map of vectors holding the pointers to the external haloed nodes.
Definition: mesh.h:138
void output_external_haloed_elements(const unsigned &p, std::ostream &outfile, const unsigned &n_plot=5)
Output all external haloed elements with processor p.
Definition: mesh.h:2204
void set_nodal_and_elemental_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with all nodal and elemental data stored in the mesh.
Definition: mesh.h:1032
virtual void setup_boundary_element_info()
Interface for function that is used to setup the boundary information (Empty virtual function – imple...
Definition: mesh.h:275
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
std::map< unsigned, Vector< GeneralisedElement * > > Root_halo_element_pt
Map of vectors holding the pointers to the root halo elements.
Definition: mesh.h:100
unsigned elemental_dimension() const
Return number of elemental dimension in mesh.
Definition: mesh.cc:8917
void set_external_haloed_node_pt(const unsigned &p, const Vector< Node * > &external_haloed_node_pt)
Set vector of external haloed node in this Mesh whose halo external counterpart is held on processor ...
Definition: mesh.h:2466
void check_inverted_elements(bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
Check for inverted elements and report outcome in boolean variable. This visits all elements at their...
Definition: mesh.cc:870
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed.
Definition: mesh.h:119
virtual void create_shared_boundaries(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, const Vector< FiniteElement * > &backed_up_f_el_pt, std::map< Data *, std::set< unsigned > > &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status)
Creates the shared boundaries, only used in unstructured meshes In this case with the "TriangleMesh" ...
Definition: mesh.h:2490
virtual void compute_norm(Vector< double > &norm)
Compute norm of solution by summing contributions of compute_norm(...) for all constituent elements i...
Definition: mesh.h:1099
void add_root_halo_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add root halo element whose non-halo counterpart is held on processor p to this Mesh.
Definition: mesh.h:1871
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
Output a given Vector function at f(n_plot) points in each element.
Definition: mesh.cc:2199
virtual void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Plot error when compared against a given exact solution. Also returns the norm of the error and that ...
Definition: mesh.h:1140
virtual void scale_mesh(const double &factor)
Scale all nodal coordinates by given factor. Virtual so it can be overloaded in SolidMesh class where...
Definition: mesh.h:373
void shift_time_values()
Shift time-dependent data along for next timestep: Deal with nodal Data/positions and the element's i...
Definition: mesh.cc:2326
Vector< GeneralisedElement * > & element_pt()
Return reference to the Vector of elements.
Definition: mesh.h:466
virtual void reorder_nodes(const bool &use_old_ordering=true)
Re-order nodes in the order in which they appear in elements – can be overloaded for more efficient r...
Definition: mesh.cc:508
void prune_halo_elements_and_nodes(Vector< GeneralisedElement * > &deleted_element_pt, const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
Definition: mesh.h:1667
unsigned nhalo_node(const unsigned &p)
Number of halo nodes in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1893
void calculate_predictions()
Calculate predictions for all Data and positions associated with the mesh, usually used in adaptive t...
Definition: mesh.cc:2366
void output_boundaries(const std::string &output_filename)
Output the nodes on the boundaries (into separate tecplot zones). Specify filename.
Definition: mesh.h:1009
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
unsigned add_external_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Add external haloed element whose non-halo counterpart is held on processor p to the storage scheme f...
Definition: mesh.cc:9475
void describe_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the dofs of the Mesh. The ostream specifies the output stream to which the descr...
Definition: mesh.cc:711
Vector< GeneralisedElement * > halo_element_pt(const unsigned &p)
Return vector of halo elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1740
void disable_output_of_halo_elements()
Function to disable halo element output.
Definition: mesh.h:2046
void output_fct_paraview(std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
Definition: mesh.cc:1491
unsigned nroot_halo_element(const unsigned &p)
Number of root halo elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1846
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:204
bool Output_halo_elements
Bool for output of halo elements.
Definition: mesh.h:2027
void build_face_mesh(const unsigned &b, Mesh *const &face_mesh_pt)
Constuct a Mesh of FACE_ELEMENTs along the b-th boundary of the mesh (which contains elements of type...
Definition: mesh.h:653
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1186
virtual void node_update(const bool &update_all_solid_nodes=false)
Update nodal positions in response to changes in the domain shape. Uses the FiniteElement::get_x(....
Definition: mesh.cc:287
void get_efficiency_of_mesh_distribution(double &av_efficiency, double &max_efficiency, double &min_efficiency)
Get efficiency of mesh distribution: In an ideal distribution without halo overhead,...
Definition: mesh.cc:6390
unsigned nroot_haloed_element()
Total number of root haloed elements in this Mesh.
Definition: mesh.h:1921
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void output_external_haloed_elements(std::ostream &outfile, const unsigned &n_plot=5)
Output all external haloed elements.
Definition: mesh.h:2190
GeneralisedElement *& root_haloed_element_pt(const unsigned &p, const unsigned &e)
Access fct to the e-th root haloed element in this Mesh whose non-halo counterpart is held on process...
Definition: mesh.h:1966
void delete_all_external_storage()
Wipe the storage for all externally-based elements.
Definition: mesh.cc:9190
unsigned nnon_halo_element()
Total number of non-halo elements in this mesh (Costly call computes result on the fly)
Definition: mesh.h:1817
virtual void compute_error(FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
Returns the norm of the error and that of the exact solution. Version with vectors of norms and error...
Definition: mesh.h:1525
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get all the halo data stored in the mesh and add pointers to the data to the map, indexed by global e...
Definition: mesh.cc:4749
Node *& external_haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th external haloed node in this Mesh whose halo external counterpart is held on p...
Definition: mesh.h:2441
unsigned nhalo_node()
Total number of halo nodes in this Mesh.
Definition: mesh.h:1878
std::map< unsigned, Vector< Node * > > Halo_node_pt
Map of vectors holding the pointers to the halo nodes.
Definition: mesh.h:106
Mesh()
Default constructor.
Definition: mesh.h:236
void set_communicator_pt(OomphCommunicator *comm_pt)
Function to set communicator (mesh is assumed to be distributed if the communicator pointer is non-nu...
Definition: mesh.h:1615
Vector< Node * > external_halo_node_pt(const unsigned &p)
Access fct to vector of external halo node in this Mesh whose non-halo external counterpart is held o...
Definition: mesh.h:2384
unsigned check_for_repeated_nodes(const double &epsilon=1.0e-12)
Check for repeated nodes within a given spatial tolerance. Return (0/1) for (pass/fail).
Definition: mesh.h:752
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this,...
Definition: mesh.cc:2590
const Vector< GeneralisedElement * > & element_pt() const
Return reference to the Vector of elements.
Definition: mesh.h:460
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
Vector< GeneralisedElement * > root_halo_element_pt(const unsigned &p)
Vector of pointers to root halo elements in this Mesh whose non-halo counterpart is held on processor...
Definition: mesh.h:1854
virtual void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Plot error when compared against a given time-depdendent exact solution. Also returns the norm of the...
Definition: mesh.h:1417
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
virtual unsigned try_to_add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Definition: mesh.h:2510
std::map< unsigned, Vector< Node * > > Haloed_node_pt
Map of vectors holding the pointers to the haloed nodes.
Definition: mesh.h:109
void doc_shared_nodes()
Doc shared nodes.
Definition: mesh.h:2072
void operator=(const Mesh &)=delete
Broken assignment operator.
void assign_initial_values_impulsive()
Assign initial values for an impulsive start.
Definition: mesh.cc:2288
Node * haloed_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th haloed node in this Mesh whose halo counterpart is held on processor p.
Definition: mesh.h:2014
bool Resize_halo_nodes_not_required
Set this to true to suppress resizing of halo nodes (at your own risk!)
Definition: mesh.h:144
void flush_node_storage()
Flush storage for nodes (only) by emptying the vectors that store the pointers to them.
Definition: mesh.h:430
void assign_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers in all elements If the boolean argument is true then also store poi...
Definition: mesh.cc:765
unsigned nexternal_halo_element(const unsigned &p)
Number of external halo elements in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:2237
Vector< Vector< Node * > > Boundary_node_pt
Vector of Vector of pointers to nodes on the boundaries: Boundary_node_pt(b,n). Note that this is pri...
Definition: mesh.h:83
virtual ~Mesh()
Virtual Destructor to clean up all memory.
Definition: mesh.cc:650
void add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Add haloed node whose halo counterpart is held on processor p to the storage scheme for haloed nodes.
Definition: mesh.h:2021
void get_external_halo_node_pt(Vector< Node * > &external_halo_node_pt)
Get vector of pointers to all external halo nodes.
Definition: mesh.h:2323
void enable_output_of_halo_elements()
Function to enable halo element output.
Definition: mesh.h:2052
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
void add_external_halo_node_pt(const unsigned &p, Node *&nod_pt)
Add external halo node whose non-halo (external) counterpart is held on processor p to the storage sc...
Definition: mesh.h:2368
void unset_keep_all_elements_as_halos()
Calll this function to unset the flag that keeps all elements in the mesh as halo elements.
Definition: mesh.h:1628
virtual void read(std::ifstream &restart_file)
Read solution from restart file.
Definition: mesh.cc:1130
void set_nodal_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Set the timestepper associated with the nodal data in the mesh.
Definition: mesh.cc:2516
void doc_mesh_distribution(DocInfo &doc_info)
Doc the mesh distribution, to be processed with tecplot macros.
Definition: mesh.cc:6470
unsigned self_test()
Self-test: Check elements and nodes. Return 0 for OK.
Definition: mesh.cc:778
void get_shared_node_pt(const unsigned &p, Vector< Node * > &shared_node_pt)
Get vector of pointers to shared nodes with processor p. Required for faster search in Missing_master...
Definition: mesh.h:2119
void get_halo_node_stats(double &av_number, unsigned &max_number, unsigned &min_number)
Get halo node stats for this distributed mesh: Average/max/min number of halo nodes over all processo...
Definition: mesh.cc:4845
void set_consistent_pinned_values_for_continuation(ContinuationStorageScheme *const &continuation_stepper_pt)
Set consistent values for pinned data in continuation.
Definition: mesh.cc:2436
Vector< Node * > prune_dead_nodes()
Prune nodes. Nodes that have been marked as obsolete are removed from the mesh (and its boundary-node...
Definition: mesh.cc:966
void output(const std::string &output_filename, const unsigned &n_plot)
Output at f(n_plot) points in each element.
Definition: mesh.h:984
std::set< int > external_halo_proc()
Return the set of processors that hold external halo nodes. This is required to avoid having to pass ...
Definition: mesh.h:2475
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
virtual void compute_error(FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Plot error when compared against a given time-dependent exact solution. Also returns the norm of the ...
Definition: mesh.h:1279
Node * get_some_non_boundary_node() const
Find a node not on any boundary in mesh_pt (useful for pinning a single node in a purely Neumann prob...
Definition: mesh.h:859
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed, i.e. if we don't have mpi).
Definition: mesh.h:1600
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
void resize_halo_nodes()
Helper function that resizes halo nodes to the same size as their non-halo counterparts if required....
Definition: mesh.cc:4426
unsigned nexternal_halo_node(const unsigned &p)
Number of external halo nodes in this Mesh whose non-halo (external) counterpart is held on processor...
Definition: mesh.h:2354
std::map< unsigned, Vector< Node * > > External_halo_node_pt
Map of vectors holding the pointers to the external halo nodes.
Definition: mesh.h:135
unsigned nhaloed_node()
Total number of haloed nodes in this Mesh.
Definition: mesh.h:1985
unsigned nexternal_haloed_node()
Total number of external haloed nodes in this Mesh.
Definition: mesh.h:2412
unsigned nshared_node(const unsigned &p)
Number of shared nodes in this Mesh who have a counterpart on processor p.
Definition: mesh.h:2097
Node * halo_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th halo node in this Mesh whose non-halo counterpart is held on processor p.
Definition: mesh.h:1914
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
Boolean used to control warning about empty mesh level timestepper function.
Definition: mesh.h:233
unsigned long assign_global_eqn_numbers(Vector< double * > &Dof_pt)
Assign the global equation numbers in the Data stored at the nodes and also internal element Data....
Definition: mesh.cc:677
virtual void compute_error(FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Plot error when compared against a given time-dependent exact solution. Also returns the norm of the ...
Definition: mesh.h:1231
virtual void distribute(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, Vector< GeneralisedElement * > &deleted_element_pt, DocInfo &doc_info, const bool &report_stats, const bool &overrule_keep_as_halo_element_status)
Distribute the problem and doc; make this virtual to allow overloading for particular meshes where fu...
Definition: mesh.cc:4959
Vector< GeneralisedElement * > haloed_element_pt(const unsigned &p)
Return vector of haloed elements in this Mesh whose haloing counterpart is held on processor p.
Definition: mesh.h:1779
Node * shared_node_pt(const unsigned &p, const unsigned &j)
Access fct to the j-th shared node in this Mesh who has a counterpart on processor p.
Definition: mesh.h:2110
void null_external_halo_node(const unsigned &p, Node *nod_pt)
Null out specified external halo node (used when deleting duplicates)
Definition: mesh.cc:8569
Mesh(const Mesh &dummy)=delete
Broken copy constructor.
unsigned nshared_node()
Total number of shared nodes in this Mesh.
Definition: mesh.h:2058
virtual void compute_error(FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Returns the norm of the error and that of the exact solution.
Definition: mesh.h:1479
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointes in external halo and haloed sch...
Definition: mesh.cc:8589
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
virtual void get_node_reordering(Vector< Node * > &reordering, const bool &use_old_ordering=true) const
Get a reordering of the nodes in the order in which they appear in elements – can be overloaded for m...
Definition: mesh.cc:532
virtual void classify_halo_and_haloed_nodes(const bool &report_stats=false)
Classify the halo and haloed nodes in the mesh. Virtual so it can be overloaded to perform additional...
Definition: mesh.h:1712
Vector< GeneralisedElement * > root_haloed_element_pt(const unsigned &p)
Vector of pointers to root haloed elements in this Mesh whose non-halo counterpart is held on process...
Definition: mesh.h:1951
void merge_meshes(const Vector< Mesh * > &sub_mesh_pt)
Merge meshes. Note: This simply merges the meshes' elements and nodes (ignoring duplicates; no bounda...
Definition: mesh.cc:65
void output_paraview(std::ofstream &file_out, const unsigned &nplot) const
Output in paraview format into specified file. Breaks up each element into sub-elements for plotting ...
Definition: mesh.cc:1225
void output_boundaries(std::ostream &outfile)
Output the nodes on the boundaries (into separate tecplot zones)
Definition: mesh.cc:1064
unsigned nexternal_halo_element()
Total number of external halo elements in this Mesh.
Definition: mesh.h:2222
virtual void dump(std::ofstream &dump_file, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
Definition: mesh.cc:1088
double total_size()
Determine the sum of all "sizes" of the FiniteElements in the mesh (non-FiniteElements are ignored)....
Definition: mesh.h:708
void dump(const std::string &dump_file_name, const bool &use_old_ordering=true) const
Dump the data in the mesh into a file for restart.
Definition: mesh.h:915
std::map< unsigned, Vector< GeneralisedElement * > > External_halo_element_pt
External halo(ed) elements are created as and when they are needed to act as source elements for the ...
Definition: mesh.h:126
Mesh(const Vector< Mesh * > &sub_mesh_pt)
Constructor builds combined mesh from the meshes specified. Note: This simply merges the meshes' elem...
Definition: mesh.h:257
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
unsigned nroot_halo_element()
Total number of root halo elements in this Mesh.
Definition: mesh.h:1830
Node * node_pt(const unsigned long &n) const
Return pointer to global node n (const version)
Definition: mesh.h:442
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
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
void resize(const unsigned &n_value)
Resize the number of equations.
Definition: nodes.cc:2167
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
General SolidMesh class.
Definition: mesh.h:2562
SolidNode * node_pt(const unsigned long &n)
Return a pointer to the n-th global SolidNode.
Definition: mesh.h:2594
void operator=(const SolidMesh &)=delete
Broken assignment operator.
SolidMesh()
Default constructor.
Definition: mesh.h:2565
void set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
Definition: mesh.cc:9564
SolidNode * boundary_node_pt(const unsigned &b, const unsigned &n)
Return n-th SolidNodes on b-th boundary.
Definition: mesh.h:2612
void scale_mesh(const double &factor)
Scale all nodal coordinates by given factor and re-assign the Lagrangian coordinates.
Definition: mesh.h:2669
SolidNode * element_node_pt(const unsigned long &e, const unsigned &n)
Return the n-th local SolidNode in elemnet e. This is required to cast the nodes in a solid mesh to b...
Definition: mesh.h:2631
static SolidICProblem Solid_IC_problem
Static problem that can be used to assign initial conditions on a given solid mesh (need to define th...
Definition: mesh.h:2679
SolidMesh(const Vector< SolidMesh * > &sub_mesh_pt)
Constructor builds combined mesh from the meshes specified. Note: This simply merges the meshes' elem...
Definition: mesh.h:2576
SolidMesh(const SolidMesh &dummy)=delete
Broken copy constructor.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:681
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
void stick_leaves_into_vector(Vector< Tree * > &)
Traverse tree and stick pointers to leaf "nodes" (only) into Vector.
Definition: tree.cc:255
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:167
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void assert_geometric_element(const unsigned &dim, const unsigned &nnode_1d=0)
Helper function to assert that finite element of type ELEMENT can be cast to base class of type GEOM_...
Definition: mesh.h:2821
bool node_global_position_comparison(Node *nd1_pt, Node *nd2_pt)
Function for ordering nodes. Return true if first node's position is "before" second nodes....
Definition: mesh.h:2918
void write_pvd_footer(std::ofstream &pvd_file)
Write the pvd file footer.
Definition: mesh.cc:9639
void write_pvd_header(std::ofstream &pvd_file)
Write the pvd file header.
Definition: mesh.cc:9615
void write_pvd_information(std::ofstream &pvd_file, const std::string &output_filename, const double &time)
Add name of output file and associated continuous time to pvd file.
Definition: mesh.cc:9624
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...