refineable_brick_mesh.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef QMESH3OOMPH_D_HEADER
#define QMESH3OOMPH_D_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include <limits.h>
 
 
// ooomph-lib includes
#include "brick_mesh.h"
#include "refineable_mesh.h"
#include "refineable_brick_element.h"
// Include to fill in additional_setup_shared_node_scheme() function
#include "refineable_mesh.template.cc"
 
namespace oomph
{
  //=======================================================================
  /// Intermediate mesh class that implements the mesh adaptation functions
  /// specified in the TreeBasedRefineableMesh class for meshes that contain the
  /// refineable variant of QElement s [The class ELEMENT provided
  /// as the template parameter must be of type
  /// RefineableQElement<3>].
  ///
  /// Mesh adaptation/refinement is implemented by OcTree
  /// procedures and any concrete implementation of this class needs to
  /// provide a OcTreeForest representation of the initial (coarse) mesh.
  //=======================================================================
  template<class ELEMENT>
  class RefineableBrickMesh : public TreeBasedRefineableMesh<ELEMENT>,
                              public virtual BrickMeshBase
  {
  public:
    /// Constructor: Setup static octree data
    RefineableBrickMesh()
    {
      // OcTree static data needs to be setup before octree-based mesh
      // refinement works
      OcTree::setup_static_data();
    }
 
    /// Broken copy constructor
    RefineableBrickMesh(const RefineableBrickMesh& dummy) = delete;
 
    /// Broken assignment operator
    void operator=(const RefineableBrickMesh&) = delete;
 
    /// Destructor:
    virtual ~RefineableBrickMesh() {}
 
    /// Set up the tree forest associated with the Mesh.
    /// Forwards call to setup_octree_forest()
    virtual void setup_tree_forest()
    {
      setup_octree_forest();
    }
 
    /// Do what it says...
    void setup_octree_forest()
    {
      if (this->Forest_pt != 0)
      {
        // Get all the tree nodes
        Vector<Tree*> all_tree_nodes_pt;
        this->Forest_pt->stick_all_tree_nodes_into_vector(all_tree_nodes_pt);
 
        // Get min and max refinement level from the tree
        unsigned local_min_ref = 0;
        unsigned local_max_ref = 0;
        this->get_refinement_levels(local_min_ref, local_max_ref);
 
        unsigned min_ref = local_min_ref;
#ifdef OOMPH_HAS_MPI
        if (Comm_pt != 0)
        {
          // Reconcile between processors: If (e.g. following
          // distribution/pruning) the mesh has no elements on this processor)
          // then ignore its contribution to the poll of max/min refinement
          // levels
          int int_local_min_ref = local_min_ref;
          if (this->nelement() == 0)
          {
            int_local_min_ref = INT_MAX;
          }
          int int_min_ref = 0;
          MPI_Allreduce(&int_local_min_ref,
                        &int_min_ref,
                        1,
                        MPI_INT,
                        MPI_MIN,
                        Comm_pt->mpi_comm());
          min_ref = int_min_ref;
        }
#endif
 
        // If we have no elements there's nothing more to be done --
        // we only came in here to participate in the communication
        if (this->nelement() == 0)
        {
          // Flush the Forest's current trees
          this->Forest_pt->flush_trees();
 
          // Delete the old Forest
          delete this->Forest_pt;
 
          // Empty dummy vector to build empty forest
          Vector<TreeRoot*> trees_pt;
 
          // Make a new (empty) Forest
          this->Forest_pt = new OcTreeForest(trees_pt);
 
          return;
        }
 
 
        // Vector to store trees for new Forest
        Vector<TreeRoot*> trees_pt;
 
        // Loop over tree nodes (e.g. elements)
        unsigned n_tree_nodes = all_tree_nodes_pt.size();
        for (unsigned e = 0; e < n_tree_nodes; e++)
        {
          Tree* tree_pt = all_tree_nodes_pt[e];
 
          // If the object_pt has been flushed then we don't want to keep
          // this tree
          if (tree_pt->object_pt() != 0)
          {
            // Get the refinement level of the current tree node
            RefineableElement* el_pt =
              dynamic_cast<RefineableElement*>(tree_pt->object_pt());
            unsigned level = el_pt->refinement_level();
 
            // If we are below the minimum refinement level, remove tree
            if (level < min_ref)
            {
              // Flush sons for this tree
              tree_pt->flush_sons();
 
              // Delete the tree (no recursion)
              delete tree_pt;
 
              // Delete the element
              delete el_pt;
            }
            else if (level == min_ref)
            {
              // Get the sons (if there are any) and store them
              unsigned n_sons = tree_pt->nsons();
              Vector<Tree*> backed_up_sons(n_sons);
              for (unsigned i_son = 0; i_son < n_sons; i_son++)
              {
                backed_up_sons[i_son] = tree_pt->son_pt(i_son);
              }
 
              // Make the element into a new treeroot
              OcTreeRoot* tree_root_pt = new OcTreeRoot(el_pt);
 
 
              // Pass sons
              tree_root_pt->set_son_pt(backed_up_sons);
 
              // Loop over sons and make the new treeroot their father
              for (unsigned i_son = 0; i_son < n_sons; i_son++)
              {
                Tree* son_pt = backed_up_sons[i_son];
 
                // Tell the son about its new father (which is also the root)
                son_pt->set_father_pt(tree_root_pt);
                son_pt->root_pt() = tree_root_pt;
 
                // ...and then tell all the descendants too
                Vector<Tree*> all_sons_pt;
                son_pt->stick_all_tree_nodes_into_vector(all_sons_pt);
                unsigned n = all_sons_pt.size();
                for (unsigned i = 0; i < n; i++)
                {
                  all_sons_pt[i]->root_pt() = tree_root_pt;
                }
              }
 
              // Add tree root to the trees_pt vector
              trees_pt.push_back(tree_root_pt);
 
              // Now kill the original (non-root) tree: First
              // flush sons for this tree
              tree_pt->flush_sons();
 
              // ...then delete the tree (no recursion)
              delete tree_pt;
            }
          }
          else // tree_pt->object_pt() is null, so delete tree
          {
            // Flush sons for this tree
            tree_pt->flush_sons();
 
            // Delete the tree (no recursion)
            delete tree_pt;
          }
        }
 
        // Flush the Forest's current trees
        this->Forest_pt->flush_trees();
 
        // Delete the old Forest
        delete this->Forest_pt;
 
        // Make a new Forest with the trees_pt roots created earlier
        this->Forest_pt = new OcTreeForest(trees_pt);
      }
      else // Create a new Forest from scratch in the "usual" uniform way
      {
        // Turn elements into individual octrees and plant in forest
        Vector<TreeRoot*> trees_pt;
        unsigned nel = nelement();
        for (unsigned iel = 0; iel < nel; iel++)
        {
          // Get pointer to full element type
          ELEMENT* el_pt = dynamic_cast<ELEMENT*>(element_pt(iel));
 
          // Build associated octree(root) -- pass pointer to corresponding
          // finite element and add the pointer to vector of octree (roots):
          OcTreeRoot* octree_root_pt = new OcTreeRoot(el_pt);
          trees_pt.push_back(octree_root_pt);
        }
        // Plant OcTreeRoots in OcTreeForest
        this->Forest_pt = new OcTreeForest(trees_pt);
      }
    }
 
  protected:
  };
 
} // namespace oomph
 
#endif