refineable_tetgen_mesh.template.cc

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//    Version 1.0; svn revision $LastChangedRevision: 1174 $
//LIC//
//LIC// $LastChangedDate: 2016-05-11 10:03:56 +0100 (Wed, 11 May 2016) $
//LIC// 
//LIC// Copyright (C) 2006-2016 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 OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC
#define OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

// OOMPH-LIB Headers
#include "refineable_tetgen_mesh.template.h"
#include "../generic/sample_point_parameters.h"
#include "../generic/mesh_as_geometric_object.h"
#include "../generic/projection.h"

namespace oomph
{

//======================================================================
/// Adapt problem based on specified elemental error estimates
//======================================================================
template <class ELEMENT>
void RefineableTetgenMesh<ELEMENT>::adapt(const Vector<double>& elem_error)
 {    


  double t_start=0.0;
  //###################################
  t_start=TimingHelpers::timer();
  //###################################

  // Get refinement targets
  Vector<double> target_size(elem_error.size());
  double max_edge_ratio=this->compute_volume_target(elem_error,
                                                    target_size);
  // Get maximum target volume
  unsigned n=target_size.size();
  double max_size=0.0;
  double min_size=DBL_MAX;
  for (unsigned e=0;e<n;e++)
   {
    if (target_size[e]>max_size) max_size=target_size[e];
    if (target_size[e]<min_size) min_size=target_size[e];
   }
  
  oomph_info << "Maximum target size: " << max_size << std::endl;
  oomph_info << "Minimum target size: " << min_size << std::endl;
  oomph_info << "Number of elements to be refined " 
             << this->Nrefined << std::endl;
  oomph_info << "Number of elements to be unrefined "
             << this->Nunrefined << std::endl;
  oomph_info << "Max edge ratio "<< max_edge_ratio << std::endl;

  double orig_max_size, orig_min_size;
  this->max_and_min_element_size(orig_max_size, orig_min_size);
  oomph_info << "Max/min element size in original mesh: " 
             << orig_max_size  << " "
             << orig_min_size << std::endl;    

  //##################################################################
  oomph_info << "adapt: Time for getting volume targets                      : "
             << TimingHelpers::timer()-t_start
             << " sec " << std::endl;
  //##################################################################

  // Should we bother to adapt?
  if ( (Nrefined > 0) || (Nunrefined > this->max_keep_unrefined()) ||
       (max_edge_ratio > this->max_permitted_edge_ratio()) )
   {

    if (! ( (Nrefined > 0) || (Nunrefined > max_keep_unrefined()) ) )
     {
      oomph_info 
       << "Mesh regeneration triggered by edge ratio criterion\n";
     }

    //###################################
    t_start=TimingHelpers::timer();
    //###################################

    // Are we dealing with a solid mesh?
    SolidMesh* solid_mesh_pt=dynamic_cast<SolidMesh*>(this);
    
    // Build temporary uniform background mesh
    //----------------------------------------
    // with volume set by maximum required volume
    //---------------------------------------
    RefineableTetgenMesh<ELEMENT>* tmp_new_mesh_pt=0;
    if (solid_mesh_pt!=0)
     {
      throw OomphLibError("Solid RefineableTetgenMesh not done yet.",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      /* tmp_new_mesh_pt=new RefineableSolidTetgenMesh<ELEMENT> */
      /*  (closed_curve_pt, */
      /*   hole_pt, */
      /*   max_size, */
      /*   this->Time_stepper_pt, */
      /*   this->Use_attributes); */
     }
    else
     {
      tmp_new_mesh_pt=new RefineableTetgenMesh<ELEMENT>
       (this->Outer_boundary_pt,
        this->Internal_surface_pt,
        max_size,
        this->Time_stepper_pt,
        this->Use_attributes,
        this->Corner_elements_must_be_split);
     }
    

    //##################################################################
    oomph_info << "adapt: Time for building temp mesh                        : "
               << TimingHelpers::timer()-t_start
               << " sec " << std::endl;
    //##################################################################


    // Get the tetgenio object associated with that mesh
    tetgenio *tmp_new_tetgenio_pt = tmp_new_mesh_pt->tetgenio_pt();
    RefineableTetgenMesh<ELEMENT>* new_mesh_pt=0;

     // If the mesh is a solid mesh then do the mapping based on the
     // Eulerian coordinates
     bool use_eulerian_coords=false;
     if (solid_mesh_pt!=0)
      {
       use_eulerian_coords=true;
      }

#ifdef OOMPH_HAS_CGAL

     // Make cgal-based bin
     CGALSamplePointContainerParameters cgal_params(this);
     if (use_eulerian_coords)
      {
       cgal_params.enable_use_eulerian_coordinates_during_setup();
      }
     MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(&cgal_params);
     
#else

     std::ostringstream error_message;
     error_message << "Non-CGAL-based target size transfer not implemented.\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);

     // Here's the relevant construction from the triangle mesh. Update...

     // // Make nonrefineable bin
     // NonRefineableBinArrayParameters params(this); 
     // if (use_eulerian_coords)
     //  {
     //   params.enable_use_eulerian_coordinates_during_setup();
     //  }
     // Vector<unsigned> bin_dim(2);
     // bin_dim[0]=Nbin_x_for_area_transfer;
     // bin_dim[1]=Nbin_y_for_area_transfer;
     // params.dimensions_of_bin_array()=bin_dim;
     // MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(&params);
     
#endif

    // Set up a map from pointer to element to its number
    // in the mesh
    std::map<GeneralisedElement*,unsigned> element_number;
    unsigned nelem=this->nelement();
    for (unsigned e=0;e<nelem;e++)
     {
      element_number[this->element_pt(e)]=e;
     }

    // Now start iterating to refine mesh recursively
    //-----------------------------------------------
    bool done=false;
    unsigned iter=0;
    while (!done)
     {
      // Accept by default but overwrite if things go wrong below
      done=true;
      
      // Loop over elements in new (tmp) mesh and visit all
      // its integration points. Check where it's located in the bin
      // structure of the current mesh and pass the target size
      // to the new element
      unsigned nelem=tmp_new_mesh_pt->nelement();
      
      // Store the target sizes for elements in the temporary
      // mesh
      Vector<double> new_transferred_target_size(nelem,0.0);
      for (unsigned e=0;e<nelem;e++)
       { 
        ELEMENT* el_pt=dynamic_cast<ELEMENT*>(tmp_new_mesh_pt->element_pt(e));
        unsigned nint=el_pt->integral_pt()->nweight();
        for (unsigned ipt=0;ipt<nint;ipt++)
         {
          // Get the coordinate of current point
          Vector<double> s(3);
          for(unsigned i=0;i<3;i++)
           {
            s[i] = el_pt->integral_pt()->knot(ipt,i);
           }
          
          Vector<double> x(3);
          el_pt->interpolated_x(s,x);
          
#if OOMPH_HAS_CGAL
          
          // Try the five nearest sample points for Newton search
          // then just settle on the nearest one
          GeomObject* geom_obj_pt=0;
          unsigned max_sample_points=5;
          dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
                                                  sample_point_container_pt())->
           limited_locate_zeta(x,max_sample_points,
                               geom_obj_pt,s);
#ifdef PARANOID
          if (geom_obj_pt==0)
           {
            std::stringstream error_message;
            error_message
             << "Limited locate zeta failed for zeta = [ "
             << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
           }
          else
           {
#endif
            FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(geom_obj_pt);
#ifdef PARANOID
            if (fe_pt==0)
             {
              std::stringstream error_message;
              error_message
               << "Cast to FE for GeomObject returned by limited locate zeta failed for zeta = [ "
               << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
              throw OomphLibError(error_message.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
             }
            else
             {
#endif
              // What's the target size of the element that contains this point
              double tg_size=target_size[element_number[fe_pt]];
              
              // Go for smallest target size over all integration 
              // points in new element
              // to force "one level" of refinement (the one-level-ness
              // is enforced below by limiting the actual reduction in
              // size
              if (new_transferred_target_size[e]!=0.0)
               {
                new_transferred_target_size[e]=
                 std::min(new_transferred_target_size[e], 
                          tg_size);
               }
              else
               {
                new_transferred_target_size[e]=tg_size;
               }
#ifdef PARANOID
             }
           }
#endif
       
// Non-CGAL   
#else
          
          std::ostringstream error_message;
          error_message 
           << "Non-CGAL-based target size transfer not implemented.\n";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);

          // Here's the relevant construction from the triangle mesh. Update...

          // // Find the bin that contains that point and its contents
          // int bin_number=0;
          // bin_array_pt->get_bin(x,bin_number);
          
          // // Did we find it?
          // if (bin_number<0)
          //  {
          //   // Not even within bin boundaries... odd
          //   std::stringstream error_message;
          //   error_message
          //    << "Very odd -- we're looking for a point[ "
          //    << x[0] << " " << x[1] << " " << x[2] << " ] that's not even \n"
          //    << "located within the bin boundaries.\n";
          //   throw OomphLibError(error_message.str(),
          //                       OOMPH_CURRENT_FUNCTION,
          //                       OOMPH_EXCEPTION_LOCATION);
          //  } // if (bin_number<0)
          // else
          //  {            
          //   // Go for smallest target size of any element in this bin
          //   // to force "one level" of refinement (the one-level-ness
          //   // is enforced below by limiting the actual reduction in
          //   // size
          //   if (new_transferred_target_size[e]!=0)
          //    {
          //      std::min(new_transferred_target_size[e],
          //               bin_min_target_size[bin_number]);
          //    }
          //   else
          //    {
          //     new_transferred_target_size[e]=bin_min_target_size[bin_number];
          //    }
            
          //  }
          
#endif
          
         } // for (ipt<nint)
        
       } // for (e<nelem)
      
      
      // do some output (keep it alive!)
      const bool output_target_sizes=false;
      if (output_target_sizes)
       {
        unsigned length=new_transferred_target_size.size();
        for (unsigned u = 0; u < length;u++)
         {
          oomph_info << "Element" << u << ",target size: "
                     << new_transferred_target_size[u] << std::endl;
         }
       }
      
      // Now copy into target size for temporary mesh but limit to
      // the equivalent of one sub-division per iteration
#ifdef OOMPH_HAS_MPI
      unsigned n_ele_need_refinement_iter = 0;
#endif 


      // Don't delete! Keep these around for debugging
      // ofstream tmp_mesh_file;
      // tmp_mesh_file.open("tmp_mesh_file.dat");
      // tmp_new_mesh_pt->output(tmp_mesh_file);
      // tmp_mesh_file.close();

      std::ofstream target_sizes_file; 
      char filename[100];
      sprintf(filename,"target_sizes%i.dat",iter);
      if (output_target_sizes)
       {
        target_sizes_file.open(filename);
       }

      const unsigned nel_new=tmp_new_mesh_pt->nelement();
      Vector<double> new_target_size(nel_new);   
      for (unsigned e=0;e<nel_new;e++)
       {
        // The finite element
        FiniteElement* f_ele_pt = tmp_new_mesh_pt->finite_element_pt(e);
        
        // Transferred target size
        const double new_size=new_transferred_target_size[e];
        if (new_size<=0.0)
         {
          std::ostringstream error_stream;
          error_stream << "This shouldn't happen! Element whose centroid is at "
                       <<  (f_ele_pt->node_pt(0)->x(0)+
                            f_ele_pt->node_pt(1)->x(0)+ 
                            f_ele_pt->node_pt(2)->x(0)+
                            f_ele_pt->node_pt(3)->x(0))/4.0 << " "
                       << (f_ele_pt->node_pt(0)->x(1)+
                           f_ele_pt->node_pt(1)->x(1)+
                           f_ele_pt->node_pt(2)->x(1)+
                           f_ele_pt->node_pt(3)->x(1))/4.0 << " "
                       << (f_ele_pt->node_pt(0)->x(2)+
                           f_ele_pt->node_pt(1)->x(2)+
                           f_ele_pt->node_pt(2)->x(2)+
                           f_ele_pt->node_pt(3)->x(2))/4.0 << " "
                       << " has no target size assigned\n";
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
         }
        else
         {
          // Limit target size to the equivalent of uniform refinement
          // during this stage of the iteration
          new_target_size[e]=new_size;
          if (new_target_size[e]<f_ele_pt->size()/4.0)
           {
            new_target_size[e]=f_ele_pt->size()/4.0;

            // ALH: It seems that tetgen "enlarges" the volume constraint
            // so this criterion can never be met unless dividing by 1.2
            // as well. MH: Is this the reason why Andrew's version of
            // adaptation never converges? Took it out.
            // new_target_size[e] /= 1.2;

            // We'll need to give it another go later
            done=false;
           } 
          
          // Don't delete! Keep around for debugging
          if (output_target_sizes)
           {
            target_sizes_file << "ZONE N=4, E=1, F=FEPOINT, ET=TETRAHEDRON\n";
            for (unsigned j=0;j<4;j++)
             {
              target_sizes_file << f_ele_pt->node_pt(j)->x(0) << " "
                                << f_ele_pt->node_pt(j)->x(1) << " "
                                << f_ele_pt->node_pt(j)->x(2) << " "
                                << new_size << " " 
                                << new_target_size[e] 
                                << std::endl;
             }
            target_sizes_file << "1 2 3 4\n"; // connectivity
            
            // Keep around; just doc at centroid
            /* target_sizes_file  */
            /*  <<  (f_ele_pt->node_pt(0)->x(0)+ */
            /*       f_ele_pt->node_pt(1)->x(0)+ */
            /*       f_ele_pt->node_pt(2)->x(0)+ */
            /*       f_ele_pt->node_pt(3)->x(0))/4.0 << " " */
            /*  << (f_ele_pt->node_pt(0)->x(1)+ */
            /*      f_ele_pt->node_pt(1)->x(1)+ */
            /*      f_ele_pt->node_pt(2)->x(1)+ */
            /*      f_ele_pt->node_pt(3)->x(1))/4.0 << " " */
            /*  << (f_ele_pt->node_pt(0)->x(2)+ */
            /*      f_ele_pt->node_pt(1)->x(2)+ */
            /*      f_ele_pt->node_pt(2)->x(2)+ */
            /*      f_ele_pt->node_pt(3)->x(2))/4.0 << " " */
            /*  << new_size << " "  */
            /*  << new_target_size[e] << std::endl; */
           }
          
#ifdef OOMPH_HAS_MPI
          // Keep track of the elements that require (un)refinement
          n_ele_need_refinement_iter++;
#endif
          
         } // else if (new_size <= 0.0)
        
       } // for (e < nel_new)
            
      // Don't delete! Keep around for debugging
      if (output_target_sizes)
       {
        target_sizes_file.close();   
       }

      if (done) 
       {
        oomph_info 
         << "All size adjustments accommodated by max. permitted size"
         << " reduction during iter " << iter << "\n";
       }
      else
       {
        oomph_info 
         << "NOT all size adjustments accommodated by max. "
         << "permitted size reduction  during iter " << iter << "\n";
       }
      

      oomph_info << "\n\n\n==================================================\n"
                 << "==================================================\n"
                 << "==================================================\n"
                 << "==================================================\n"
                 << "\n\n\n";

      //##################################################################
      oomph_info << "adapt: Time for new_target_size[.]                      : "
                 << TimingHelpers::timer()-t_start
                 << " sec " << std::endl;
      //##################################################################


      // Now create the new mesh from TriangulateIO structure
      //-----------------------------------------------------
      // associated with uniform background mesh and the
      //------------------------------------------------
      // associated target element sizes.
      //---------------------------------
      
      //###################################
      t_start=TimingHelpers::timer();
      //###################################
      
      // Solid mesh?
      if (solid_mesh_pt!=0)
       {
        std::ostringstream error_message;
        error_message 
         << "RefineableSolidTetgenMesh not implemented yet.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
        /* new_mesh_pt=new RefineableSolidTetgenMesh<ELEMENT> */
        /*  (new_target_area, */
        /*   tmp_new_triangulateio, */
        /*   this->Time_stepper_pt, */
        /*   this->Use_attributes); */
       }      
      // No solid mesh
      else
       { 
        new_mesh_pt=new RefineableTetgenMesh<ELEMENT>
         (new_target_size,
          tmp_new_tetgenio_pt,
          this->Outer_boundary_pt,
          this->Internal_surface_pt,
          this->Time_stepper_pt,
          this->Use_attributes);
       }    
      
      //##################################################################
      oomph_info << "adapt: Time for new_mesh_pt                            : "
                 << TimingHelpers::timer()-t_start
                 << " sec " << std::endl;
      //##################################################################
      
      
      // Not done: get ready for another iteration
      iter++;


      // Check if the new mesh actually differs from the old one
      // If not, we're done.
      if (!done)
       {
        unsigned nnod=tmp_new_mesh_pt->nnode();
        if (nnod==new_mesh_pt->nnode())
         {
          unsigned nel=tmp_new_mesh_pt->nelement();
          if (nel==new_mesh_pt->nelement())
           {
            bool nodes_identical=true;
            for (unsigned j=0;j<nnod;j++)
             {
              bool coords_identical=true;
              for (unsigned i=0;i<3;i++)
               {
                if (new_mesh_pt->node_pt(j)->x(i)!=
                    tmp_new_mesh_pt->node_pt(j)->x(i))
                 {
                  coords_identical=false;
                 }
               }
              if (!coords_identical)
               {
                nodes_identical=false;                    
                break;
               }
             }
            if (nodes_identical)
             {
              done=true;
             }
           }
         }
       }

      //Delete the temporary mesh
      delete tmp_new_mesh_pt;

      //Now transfer over the pointers
      if (!done)
       {
        tmp_new_mesh_pt=new_mesh_pt;
        tmp_new_tetgenio_pt=new_mesh_pt->tetgenio_pt();
       }
      
     } // end of iteration
    

    // Check that the projection step is not disabled
    if (!Projection_is_disabled)
     {
      //###################################
      t_start=TimingHelpers::timer();
      //###################################
      
      // Project current solution onto new mesh
      //---------------------------------------
      ProjectionProblem<ELEMENT>* project_problem_pt=
       new ProjectionProblem<ELEMENT>;
      project_problem_pt->mesh_pt()=new_mesh_pt;
      project_problem_pt->project(this);
      delete project_problem_pt;

      //##################################################################
      oomph_info 
       << "adapt: Time for project soln onto new mesh                : "
       << TimingHelpers::timer()-t_start
       << " sec " << std::endl;
      //##################################################################
     }

    //###################################
    t_start=TimingHelpers::timer();
    //###################################
    
    //this->output("pre_proj",5);
    //new_mesh_pt->output("post_proj.dat",5);
    
    //Flush the old mesh 
    unsigned nnod=nnode();
    for(unsigned j=nnod;j>0;j--)  
     { 
      delete Node_pt[j-1];  
      Node_pt[j-1] = 0; 
     } 
    unsigned nel=nelement(); 
    for(unsigned e=nel;e>0;e--)  
     { 
      delete Element_pt[e-1];  
      Element_pt[e-1] = 0; 
     } 
    
    // Now copy back to current mesh
    //------------------------------
    nnod=new_mesh_pt->nnode();
    Node_pt.resize(nnod);
    nel=new_mesh_pt->nelement();
    Element_pt.resize(nel);  
    for(unsigned j=0;j<nnod;j++)
     { 
      Node_pt[j] = new_mesh_pt->node_pt(j);
     } 
    for(unsigned e=0;e<nel;e++)
     { 
      Element_pt[e] = new_mesh_pt->element_pt(e);
     } 
    
    //Copy the boundary schemes
    unsigned nbound=new_mesh_pt->nboundary();
    Boundary_element_pt.resize(nbound);
    Face_index_at_boundary.resize(nbound);
    Boundary_node_pt.resize(nbound);
    for (unsigned b=0;b<nbound;b++)
     {
      unsigned nel=new_mesh_pt->nboundary_element(b);
      Boundary_element_pt[b].resize(nel);
      Face_index_at_boundary[b].resize(nel);
      for (unsigned e=0;e<nel;e++)
       {
        Boundary_element_pt[b][e]=new_mesh_pt->boundary_element_pt(b,e);
        Face_index_at_boundary[b][e]=new_mesh_pt->face_index_at_boundary(b,e);
       }
      unsigned nnod=new_mesh_pt->nboundary_node(b);
      Boundary_node_pt[b].resize(nnod);
      for (unsigned j=0;j<nnod;j++)
       {
        Boundary_node_pt[b][j]=new_mesh_pt->boundary_node_pt(b,j);
       }
     }

    //Also copy over the new boundary and region information
    unsigned n_region = new_mesh_pt->nregion();

    //Only bother if we have regions
    if(n_region > 1)
     {
      //Deal with the region information first
      this->Region_element_pt.resize(n_region);
      this->Region_attribute.resize(n_region);
      for(unsigned i=0;i<n_region;i++)
       {
        // Copy across region attributes (region ids!)
        this->Region_attribute[i] = new_mesh_pt->region_attribute(i);

        //Find the number of elements in the region
        unsigned r=this->Region_attribute[i];
        unsigned n_region_element = new_mesh_pt->nregion_element(r);
        this->Region_element_pt[i].resize(n_region_element);
        for(unsigned e=0;e<n_region_element;e++)
         {
          this->Region_element_pt[i][e] = new_mesh_pt->region_element_pt(r,e);
         }
       }

      //Now the boundary region information
      this->Boundary_region_element_pt.resize(nbound);
      this->Face_index_region_at_boundary.resize(nbound);
      
      //Now loop over the boundaries
      for(unsigned b=0;b<nbound;++b)
       {
        //Loop over the regions
        for(unsigned i=0;i<n_region;++i)
         {        
          unsigned r=this->Region_attribute[i];

          unsigned n_boundary_el_in_region = 
           new_mesh_pt->nboundary_element_in_region(b,r);
          if(n_boundary_el_in_region > 0)
           {
            //Allocate storage in the map
            this->Boundary_region_element_pt[b][r].
             resize(n_boundary_el_in_region);
            this->Face_index_region_at_boundary[b][r].
             resize(n_boundary_el_in_region);

            //Copy over the information
            for(unsigned e=0;e<n_boundary_el_in_region;++e)
             {
              this->Boundary_region_element_pt[b][r][e]
               = new_mesh_pt->boundary_element_in_region_pt(b,r,e);
              this->Face_index_region_at_boundary[b][r][e] 
               = new_mesh_pt->face_index_at_boundary_in_region(b,r,e);
             }
           }
         }
       } //End of loop over boundaries

     } //End of case when more than one region
    
    // Copy TriangulateIO representation
    this->set_deep_copy_tetgenio_pt(new_mesh_pt->tetgenio_pt());
    
    // Flush the mesh
    new_mesh_pt->flush_element_and_node_storage();
    
    // Delete the mesh and the problem
    delete new_mesh_pt;


    //##################################################################
    oomph_info << "adapt: Time for moving nodes etc. to actual mesh          : "
               << TimingHelpers::timer()-t_start
               << " sec " << std::endl;
    //##################################################################

    // Solid mesh?
    if (solid_mesh_pt!=0)
     {
      // Warning
      std::stringstream error_message;
      error_message 
       << "Lagrangian coordinates are currently not projected but are\n"
       << "are re-set during adaptation. This is not appropriate for\n"
       << "real solid mechanics problems!\n";
      OomphLibWarning(error_message.str(),
                      OOMPH_CURRENT_FUNCTION,
                      OOMPH_EXCEPTION_LOCATION);
      
      // Reset Lagrangian coordinates
      dynamic_cast<SolidMesh*>(this)->set_lagrangian_nodal_coordinates();
     }
    
    double max_area;
    double min_area;
    this->max_and_min_element_size(max_area, min_area);
    oomph_info << "Max/min element size in adapted mesh: " 
               << max_area  << " "
               << min_area << std::endl;    
   }
  else
   {
    oomph_info << "Not enough benefit in adaptation.\n";
    Nrefined=0;
    Nunrefined=0;
   }
 }

}

#endif