refineable_tetgen_mesh.template.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1174 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-11 10:03:56 +0100 (Wed, 11 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #ifndef OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC
32 #define OOMPH_REFINEABLE_TETGEN_MESH_TEMPLATE_CC
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 // OOMPH-LIB Headers
40 #include "refineable_tetgen_mesh.template.h"
41 #include "../generic/sample_point_parameters.h"
42 #include "../generic/mesh_as_geometric_object.h"
43 #include "../generic/projection.h"
44 
45 namespace oomph
46 {
47 
48 //======================================================================
49 /// Adapt problem based on specified elemental error estimates
50 //======================================================================
51 template <class ELEMENT>
52 void RefineableTetgenMesh<ELEMENT>::adapt(const Vector<double>& elem_error)
53  {
54 
55 
56  double t_start=0.0;
57  //###################################
58  t_start=TimingHelpers::timer();
59  //###################################
60 
61  // Get refinement targets
62  Vector<double> target_size(elem_error.size());
63  double max_edge_ratio=this->compute_volume_target(elem_error,
64  target_size);
65  // Get maximum target volume
66  unsigned n=target_size.size();
67  double max_size=0.0;
68  double min_size=DBL_MAX;
69  for (unsigned e=0;e<n;e++)
70  {
71  if (target_size[e]>max_size) max_size=target_size[e];
72  if (target_size[e]<min_size) min_size=target_size[e];
73  }
74 
75  oomph_info << "Maximum target size: " << max_size << std::endl;
76  oomph_info << "Minimum target size: " << min_size << std::endl;
77  oomph_info << "Number of elements to be refined "
78  << this->Nrefined << std::endl;
79  oomph_info << "Number of elements to be unrefined "
80  << this->Nunrefined << std::endl;
81  oomph_info << "Max edge ratio "<< max_edge_ratio << std::endl;
82 
83  double orig_max_size, orig_min_size;
84  this->max_and_min_element_size(orig_max_size, orig_min_size);
85  oomph_info << "Max/min element size in original mesh: "
86  << orig_max_size << " "
87  << orig_min_size << std::endl;
88 
89  //##################################################################
90  oomph_info << "adapt: Time for getting volume targets : "
91  << TimingHelpers::timer()-t_start
92  << " sec " << std::endl;
93  //##################################################################
94 
95  // Should we bother to adapt?
96  if ( (Nrefined > 0) || (Nunrefined > this->max_keep_unrefined()) ||
97  (max_edge_ratio > this->max_permitted_edge_ratio()) )
98  {
99 
100  if (! ( (Nrefined > 0) || (Nunrefined > max_keep_unrefined()) ) )
101  {
102  oomph_info
103  << "Mesh regeneration triggered by edge ratio criterion\n";
104  }
105 
106  //###################################
107  t_start=TimingHelpers::timer();
108  //###################################
109 
110  // Are we dealing with a solid mesh?
111  SolidMesh* solid_mesh_pt=dynamic_cast<SolidMesh*>(this);
112 
113  // Build temporary uniform background mesh
114  //----------------------------------------
115  // with volume set by maximum required volume
116  //---------------------------------------
117  RefineableTetgenMesh<ELEMENT>* tmp_new_mesh_pt=0;
118  if (solid_mesh_pt!=0)
119  {
120  throw OomphLibError("Solid RefineableTetgenMesh not done yet.",
121  OOMPH_CURRENT_FUNCTION,
122  OOMPH_EXCEPTION_LOCATION);
123  /* tmp_new_mesh_pt=new RefineableSolidTetgenMesh<ELEMENT> */
124  /* (closed_curve_pt, */
125  /* hole_pt, */
126  /* max_size, */
127  /* this->Time_stepper_pt, */
128  /* this->Use_attributes); */
129  }
130  else
131  {
132  tmp_new_mesh_pt=new RefineableTetgenMesh<ELEMENT>
133  (this->Outer_boundary_pt,
134  this->Internal_surface_pt,
135  max_size,
136  this->Time_stepper_pt,
137  this->Use_attributes,
138  this->Corner_elements_must_be_split);
139  }
140 
141 
142  //##################################################################
143  oomph_info << "adapt: Time for building temp mesh : "
144  << TimingHelpers::timer()-t_start
145  << " sec " << std::endl;
146  //##################################################################
147 
148 
149  // Get the tetgenio object associated with that mesh
150  tetgenio *tmp_new_tetgenio_pt = tmp_new_mesh_pt->tetgenio_pt();
151  RefineableTetgenMesh<ELEMENT>* new_mesh_pt=0;
152 
153  // If the mesh is a solid mesh then do the mapping based on the
154  // Eulerian coordinates
155  bool use_eulerian_coords=false;
156  if (solid_mesh_pt!=0)
157  {
158  use_eulerian_coords=true;
159  }
160 
161 #ifdef OOMPH_HAS_CGAL
162 
163  // Make cgal-based bin
164  CGALSamplePointContainerParameters cgal_params(this);
165  if (use_eulerian_coords)
166  {
167  cgal_params.enable_use_eulerian_coordinates_during_setup();
168  }
169  MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(&cgal_params);
170 
171 #else
172 
173  std::ostringstream error_message;
174  error_message << "Non-CGAL-based target size transfer not implemented.\n";
175  throw OomphLibError(error_message.str(),
176  OOMPH_CURRENT_FUNCTION,
177  OOMPH_EXCEPTION_LOCATION);
178 
179  // Here's the relevant construction from the triangle mesh. Update...
180 
181  // // Make nonrefineable bin
182  // NonRefineableBinArrayParameters params(this);
183  // if (use_eulerian_coords)
184  // {
185  // params.enable_use_eulerian_coordinates_during_setup();
186  // }
187  // Vector<unsigned> bin_dim(2);
188  // bin_dim[0]=Nbin_x_for_area_transfer;
189  // bin_dim[1]=Nbin_y_for_area_transfer;
190  // params.dimensions_of_bin_array()=bin_dim;
191  // MeshAsGeomObject* mesh_geom_obj_pt=new MeshAsGeomObject(&params);
192 
193 #endif
194 
195  // Set up a map from pointer to element to its number
196  // in the mesh
197  std::map<GeneralisedElement*,unsigned> element_number;
198  unsigned nelem=this->nelement();
199  for (unsigned e=0;e<nelem;e++)
200  {
201  element_number[this->element_pt(e)]=e;
202  }
203 
204  // Now start iterating to refine mesh recursively
205  //-----------------------------------------------
206  bool done=false;
207  unsigned iter=0;
208  while (!done)
209  {
210  // Accept by default but overwrite if things go wrong below
211  done=true;
212 
213  // Loop over elements in new (tmp) mesh and visit all
214  // its integration points. Check where it's located in the bin
215  // structure of the current mesh and pass the target size
216  // to the new element
217  unsigned nelem=tmp_new_mesh_pt->nelement();
218 
219  // Store the target sizes for elements in the temporary
220  // mesh
221  Vector<double> new_transferred_target_size(nelem,0.0);
222  for (unsigned e=0;e<nelem;e++)
223  {
224  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(tmp_new_mesh_pt->element_pt(e));
225  unsigned nint=el_pt->integral_pt()->nweight();
226  for (unsigned ipt=0;ipt<nint;ipt++)
227  {
228  // Get the coordinate of current point
229  Vector<double> s(3);
230  for(unsigned i=0;i<3;i++)
231  {
232  s[i] = el_pt->integral_pt()->knot(ipt,i);
233  }
234 
235  Vector<double> x(3);
236  el_pt->interpolated_x(s,x);
237 
238 #if OOMPH_HAS_CGAL
239 
240  // Try the five nearest sample points for Newton search
241  // then just settle on the nearest one
242  GeomObject* geom_obj_pt=0;
243  unsigned max_sample_points=5;
244  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
245  sample_point_container_pt())->
246  limited_locate_zeta(x,max_sample_points,
247  geom_obj_pt,s);
248 #ifdef PARANOID
249  if (geom_obj_pt==0)
250  {
251  std::stringstream error_message;
252  error_message
253  << "Limited locate zeta failed for zeta = [ "
254  << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
255  throw OomphLibError(error_message.str(),
256  OOMPH_CURRENT_FUNCTION,
257  OOMPH_EXCEPTION_LOCATION);
258  }
259  else
260  {
261 #endif
262  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(geom_obj_pt);
263 #ifdef PARANOID
264  if (fe_pt==0)
265  {
266  std::stringstream error_message;
267  error_message
268  << "Cast to FE for GeomObject returned by limited locate zeta failed for zeta = [ "
269  << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
270  throw OomphLibError(error_message.str(),
271  OOMPH_CURRENT_FUNCTION,
272  OOMPH_EXCEPTION_LOCATION);
273  }
274  else
275  {
276 #endif
277  // What's the target size of the element that contains this point
278  double tg_size=target_size[element_number[fe_pt]];
279 
280  // Go for smallest target size over all integration
281  // points in new element
282  // to force "one level" of refinement (the one-level-ness
283  // is enforced below by limiting the actual reduction in
284  // size
285  if (new_transferred_target_size[e]!=0.0)
286  {
287  new_transferred_target_size[e]=
288  std::min(new_transferred_target_size[e],
289  tg_size);
290  }
291  else
292  {
293  new_transferred_target_size[e]=tg_size;
294  }
295 #ifdef PARANOID
296  }
297  }
298 #endif
299 
300 // Non-CGAL
301 #else
302 
303  std::ostringstream error_message;
304  error_message
305  << "Non-CGAL-based target size transfer not implemented.\n";
306  throw OomphLibError(error_message.str(),
307  OOMPH_CURRENT_FUNCTION,
308  OOMPH_EXCEPTION_LOCATION);
309 
310  // Here's the relevant construction from the triangle mesh. Update...
311 
312  // // Find the bin that contains that point and its contents
313  // int bin_number=0;
314  // bin_array_pt->get_bin(x,bin_number);
315 
316  // // Did we find it?
317  // if (bin_number<0)
318  // {
319  // // Not even within bin boundaries... odd
320  // std::stringstream error_message;
321  // error_message
322  // << "Very odd -- we're looking for a point[ "
323  // << x[0] << " " << x[1] << " " << x[2] << " ] that's not even \n"
324  // << "located within the bin boundaries.\n";
325  // throw OomphLibError(error_message.str(),
326  // OOMPH_CURRENT_FUNCTION,
327  // OOMPH_EXCEPTION_LOCATION);
328  // } // if (bin_number<0)
329  // else
330  // {
331  // // Go for smallest target size of any element in this bin
332  // // to force "one level" of refinement (the one-level-ness
333  // // is enforced below by limiting the actual reduction in
334  // // size
335  // if (new_transferred_target_size[e]!=0)
336  // {
337  // std::min(new_transferred_target_size[e],
338  // bin_min_target_size[bin_number]);
339  // }
340  // else
341  // {
342  // new_transferred_target_size[e]=bin_min_target_size[bin_number];
343  // }
344 
345  // }
346 
347 #endif
348 
349  } // for (ipt<nint)
350 
351  } // for (e<nelem)
352 
353 
354  // do some output (keep it alive!)
355  const bool output_target_sizes=false;
356  if (output_target_sizes)
357  {
358  unsigned length=new_transferred_target_size.size();
359  for (unsigned u = 0; u < length;u++)
360  {
361  oomph_info << "Element" << u << ",target size: "
362  << new_transferred_target_size[u] << std::endl;
363  }
364  }
365 
366  // Now copy into target size for temporary mesh but limit to
367  // the equivalent of one sub-division per iteration
368 #ifdef OOMPH_HAS_MPI
369  unsigned n_ele_need_refinement_iter = 0;
370 #endif
371 
372 
373  // Don't delete! Keep these around for debugging
374  // ofstream tmp_mesh_file;
375  // tmp_mesh_file.open("tmp_mesh_file.dat");
376  // tmp_new_mesh_pt->output(tmp_mesh_file);
377  // tmp_mesh_file.close();
378 
379  std::ofstream target_sizes_file;
380  char filename[100];
381  sprintf(filename,"target_sizes%i.dat",iter);
382  if (output_target_sizes)
383  {
384  target_sizes_file.open(filename);
385  }
386 
387  const unsigned nel_new=tmp_new_mesh_pt->nelement();
388  Vector<double> new_target_size(nel_new);
389  for (unsigned e=0;e<nel_new;e++)
390  {
391  // The finite element
392  FiniteElement* f_ele_pt = tmp_new_mesh_pt->finite_element_pt(e);
393 
394  // Transferred target size
395  const double new_size=new_transferred_target_size[e];
396  if (new_size<=0.0)
397  {
398  std::ostringstream error_stream;
399  error_stream << "This shouldn't happen! Element whose centroid is at "
400  << (f_ele_pt->node_pt(0)->x(0)+
401  f_ele_pt->node_pt(1)->x(0)+
402  f_ele_pt->node_pt(2)->x(0)+
403  f_ele_pt->node_pt(3)->x(0))/4.0 << " "
404  << (f_ele_pt->node_pt(0)->x(1)+
405  f_ele_pt->node_pt(1)->x(1)+
406  f_ele_pt->node_pt(2)->x(1)+
407  f_ele_pt->node_pt(3)->x(1))/4.0 << " "
408  << (f_ele_pt->node_pt(0)->x(2)+
409  f_ele_pt->node_pt(1)->x(2)+
410  f_ele_pt->node_pt(2)->x(2)+
411  f_ele_pt->node_pt(3)->x(2))/4.0 << " "
412  << " has no target size assigned\n";
413  throw OomphLibError(error_stream.str(),
414  OOMPH_CURRENT_FUNCTION,
415  OOMPH_EXCEPTION_LOCATION);
416  }
417  else
418  {
419  // Limit target size to the equivalent of uniform refinement
420  // during this stage of the iteration
421  new_target_size[e]=new_size;
422  if (new_target_size[e]<f_ele_pt->size()/4.0)
423  {
424  new_target_size[e]=f_ele_pt->size()/4.0;
425 
426  // ALH: It seems that tetgen "enlarges" the volume constraint
427  // so this criterion can never be met unless dividing by 1.2
428  // as well. MH: Is this the reason why Andrew's version of
429  // adaptation never converges? Took it out.
430  // new_target_size[e] /= 1.2;
431 
432  // We'll need to give it another go later
433  done=false;
434  }
435 
436  // Don't delete! Keep around for debugging
437  if (output_target_sizes)
438  {
439  target_sizes_file << "ZONE N=4, E=1, F=FEPOINT, ET=TETRAHEDRON\n";
440  for (unsigned j=0;j<4;j++)
441  {
442  target_sizes_file << f_ele_pt->node_pt(j)->x(0) << " "
443  << f_ele_pt->node_pt(j)->x(1) << " "
444  << f_ele_pt->node_pt(j)->x(2) << " "
445  << new_size << " "
446  << new_target_size[e]
447  << std::endl;
448  }
449  target_sizes_file << "1 2 3 4\n"; // connectivity
450 
451  // Keep around; just doc at centroid
452  /* target_sizes_file */
453  /* << (f_ele_pt->node_pt(0)->x(0)+ */
454  /* f_ele_pt->node_pt(1)->x(0)+ */
455  /* f_ele_pt->node_pt(2)->x(0)+ */
456  /* f_ele_pt->node_pt(3)->x(0))/4.0 << " " */
457  /* << (f_ele_pt->node_pt(0)->x(1)+ */
458  /* f_ele_pt->node_pt(1)->x(1)+ */
459  /* f_ele_pt->node_pt(2)->x(1)+ */
460  /* f_ele_pt->node_pt(3)->x(1))/4.0 << " " */
461  /* << (f_ele_pt->node_pt(0)->x(2)+ */
462  /* f_ele_pt->node_pt(1)->x(2)+ */
463  /* f_ele_pt->node_pt(2)->x(2)+ */
464  /* f_ele_pt->node_pt(3)->x(2))/4.0 << " " */
465  /* << new_size << " " */
466  /* << new_target_size[e] << std::endl; */
467  }
468 
469 #ifdef OOMPH_HAS_MPI
470  // Keep track of the elements that require (un)refinement
471  n_ele_need_refinement_iter++;
472 #endif
473 
474  } // else if (new_size <= 0.0)
475 
476  } // for (e < nel_new)
477 
478  // Don't delete! Keep around for debugging
479  if (output_target_sizes)
480  {
481  target_sizes_file.close();
482  }
483 
484  if (done)
485  {
486  oomph_info
487  << "All size adjustments accommodated by max. permitted size"
488  << " reduction during iter " << iter << "\n";
489  }
490  else
491  {
492  oomph_info
493  << "NOT all size adjustments accommodated by max. "
494  << "permitted size reduction during iter " << iter << "\n";
495  }
496 
497 
498  oomph_info << "\n\n\n==================================================\n"
499  << "==================================================\n"
500  << "==================================================\n"
501  << "==================================================\n"
502  << "\n\n\n";
503 
504  //##################################################################
505  oomph_info << "adapt: Time for new_target_size[.] : "
506  << TimingHelpers::timer()-t_start
507  << " sec " << std::endl;
508  //##################################################################
509 
510 
511  // Now create the new mesh from TriangulateIO structure
512  //-----------------------------------------------------
513  // associated with uniform background mesh and the
514  //------------------------------------------------
515  // associated target element sizes.
516  //---------------------------------
517 
518  //###################################
519  t_start=TimingHelpers::timer();
520  //###################################
521 
522  // Solid mesh?
523  if (solid_mesh_pt!=0)
524  {
525  std::ostringstream error_message;
526  error_message
527  << "RefineableSolidTetgenMesh not implemented yet.\n";
528  throw OomphLibError(error_message.str(),
529  OOMPH_CURRENT_FUNCTION,
530  OOMPH_EXCEPTION_LOCATION);
531  /* new_mesh_pt=new RefineableSolidTetgenMesh<ELEMENT> */
532  /* (new_target_area, */
533  /* tmp_new_triangulateio, */
534  /* this->Time_stepper_pt, */
535  /* this->Use_attributes); */
536  }
537  // No solid mesh
538  else
539  {
540  new_mesh_pt=new RefineableTetgenMesh<ELEMENT>
541  (new_target_size,
542  tmp_new_tetgenio_pt,
543  this->Outer_boundary_pt,
544  this->Internal_surface_pt,
545  this->Time_stepper_pt,
546  this->Use_attributes);
547  }
548 
549  //##################################################################
550  oomph_info << "adapt: Time for new_mesh_pt : "
551  << TimingHelpers::timer()-t_start
552  << " sec " << std::endl;
553  //##################################################################
554 
555 
556  // Not done: get ready for another iteration
557  iter++;
558 
559 
560  // Check if the new mesh actually differs from the old one
561  // If not, we're done.
562  if (!done)
563  {
564  unsigned nnod=tmp_new_mesh_pt->nnode();
565  if (nnod==new_mesh_pt->nnode())
566  {
567  unsigned nel=tmp_new_mesh_pt->nelement();
568  if (nel==new_mesh_pt->nelement())
569  {
570  bool nodes_identical=true;
571  for (unsigned j=0;j<nnod;j++)
572  {
573  bool coords_identical=true;
574  for (unsigned i=0;i<3;i++)
575  {
576  if (new_mesh_pt->node_pt(j)->x(i)!=
577  tmp_new_mesh_pt->node_pt(j)->x(i))
578  {
579  coords_identical=false;
580  }
581  }
582  if (!coords_identical)
583  {
584  nodes_identical=false;
585  break;
586  }
587  }
588  if (nodes_identical)
589  {
590  done=true;
591  }
592  }
593  }
594  }
595 
596  //Delete the temporary mesh
597  delete tmp_new_mesh_pt;
598 
599  //Now transfer over the pointers
600  if (!done)
601  {
602  tmp_new_mesh_pt=new_mesh_pt;
603  tmp_new_tetgenio_pt=new_mesh_pt->tetgenio_pt();
604  }
605 
606  } // end of iteration
607 
608 
609  // Check that the projection step is not disabled
610  if (!Projection_is_disabled)
611  {
612  //###################################
613  t_start=TimingHelpers::timer();
614  //###################################
615 
616  // Project current solution onto new mesh
617  //---------------------------------------
618  ProjectionProblem<ELEMENT>* project_problem_pt=
619  new ProjectionProblem<ELEMENT>;
620  project_problem_pt->mesh_pt()=new_mesh_pt;
621  project_problem_pt->project(this);
622  delete project_problem_pt;
623 
624  //##################################################################
625  oomph_info
626  << "adapt: Time for project soln onto new mesh : "
627  << TimingHelpers::timer()-t_start
628  << " sec " << std::endl;
629  //##################################################################
630  }
631 
632  //###################################
633  t_start=TimingHelpers::timer();
634  //###################################
635 
636  //this->output("pre_proj",5);
637  //new_mesh_pt->output("post_proj.dat",5);
638 
639  //Flush the old mesh
640  unsigned nnod=nnode();
641  for(unsigned j=nnod;j>0;j--)
642  {
643  delete Node_pt[j-1];
644  Node_pt[j-1] = 0;
645  }
646  unsigned nel=nelement();
647  for(unsigned e=nel;e>0;e--)
648  {
649  delete Element_pt[e-1];
650  Element_pt[e-1] = 0;
651  }
652 
653  // Now copy back to current mesh
654  //------------------------------
655  nnod=new_mesh_pt->nnode();
656  Node_pt.resize(nnod);
657  nel=new_mesh_pt->nelement();
658  Element_pt.resize(nel);
659  for(unsigned j=0;j<nnod;j++)
660  {
661  Node_pt[j] = new_mesh_pt->node_pt(j);
662  }
663  for(unsigned e=0;e<nel;e++)
664  {
665  Element_pt[e] = new_mesh_pt->element_pt(e);
666  }
667 
668  //Copy the boundary schemes
669  unsigned nbound=new_mesh_pt->nboundary();
670  Boundary_element_pt.resize(nbound);
671  Face_index_at_boundary.resize(nbound);
672  Boundary_node_pt.resize(nbound);
673  for (unsigned b=0;b<nbound;b++)
674  {
675  unsigned nel=new_mesh_pt->nboundary_element(b);
676  Boundary_element_pt[b].resize(nel);
677  Face_index_at_boundary[b].resize(nel);
678  for (unsigned e=0;e<nel;e++)
679  {
680  Boundary_element_pt[b][e]=new_mesh_pt->boundary_element_pt(b,e);
681  Face_index_at_boundary[b][e]=new_mesh_pt->face_index_at_boundary(b,e);
682  }
683  unsigned nnod=new_mesh_pt->nboundary_node(b);
684  Boundary_node_pt[b].resize(nnod);
685  for (unsigned j=0;j<nnod;j++)
686  {
687  Boundary_node_pt[b][j]=new_mesh_pt->boundary_node_pt(b,j);
688  }
689  }
690 
691  //Also copy over the new boundary and region information
692  unsigned n_region = new_mesh_pt->nregion();
693 
694  //Only bother if we have regions
695  if(n_region > 1)
696  {
697  //Deal with the region information first
698  this->Region_element_pt.resize(n_region);
699  this->Region_attribute.resize(n_region);
700  for(unsigned i=0;i<n_region;i++)
701  {
702  // Copy across region attributes (region ids!)
703  this->Region_attribute[i] = new_mesh_pt->region_attribute(i);
704 
705  //Find the number of elements in the region
706  unsigned r=this->Region_attribute[i];
707  unsigned n_region_element = new_mesh_pt->nregion_element(r);
708  this->Region_element_pt[i].resize(n_region_element);
709  for(unsigned e=0;e<n_region_element;e++)
710  {
711  this->Region_element_pt[i][e] = new_mesh_pt->region_element_pt(r,e);
712  }
713  }
714 
715  //Now the boundary region information
716  this->Boundary_region_element_pt.resize(nbound);
717  this->Face_index_region_at_boundary.resize(nbound);
718 
719  //Now loop over the boundaries
720  for(unsigned b=0;b<nbound;++b)
721  {
722  //Loop over the regions
723  for(unsigned i=0;i<n_region;++i)
724  {
725  unsigned r=this->Region_attribute[i];
726 
727  unsigned n_boundary_el_in_region =
728  new_mesh_pt->nboundary_element_in_region(b,r);
729  if(n_boundary_el_in_region > 0)
730  {
731  //Allocate storage in the map
732  this->Boundary_region_element_pt[b][r].
733  resize(n_boundary_el_in_region);
734  this->Face_index_region_at_boundary[b][r].
735  resize(n_boundary_el_in_region);
736 
737  //Copy over the information
738  for(unsigned e=0;e<n_boundary_el_in_region;++e)
739  {
740  this->Boundary_region_element_pt[b][r][e]
741  = new_mesh_pt->boundary_element_in_region_pt(b,r,e);
742  this->Face_index_region_at_boundary[b][r][e]
743  = new_mesh_pt->face_index_at_boundary_in_region(b,r,e);
744  }
745  }
746  }
747  } //End of loop over boundaries
748 
749  } //End of case when more than one region
750 
751  // Copy TriangulateIO representation
752  this->set_deep_copy_tetgenio_pt(new_mesh_pt->tetgenio_pt());
753 
754  // Flush the mesh
755  new_mesh_pt->flush_element_and_node_storage();
756 
757  // Delete the mesh and the problem
758  delete new_mesh_pt;
759 
760 
761  //##################################################################
762  oomph_info << "adapt: Time for moving nodes etc. to actual mesh : "
763  << TimingHelpers::timer()-t_start
764  << " sec " << std::endl;
765  //##################################################################
766 
767  // Solid mesh?
768  if (solid_mesh_pt!=0)
769  {
770  // Warning
771  std::stringstream error_message;
772  error_message
773  << "Lagrangian coordinates are currently not projected but are\n"
774  << "are re-set during adaptation. This is not appropriate for\n"
775  << "real solid mechanics problems!\n";
776  OomphLibWarning(error_message.str(),
777  OOMPH_CURRENT_FUNCTION,
778  OOMPH_EXCEPTION_LOCATION);
779 
780  // Reset Lagrangian coordinates
781  dynamic_cast<SolidMesh*>(this)->set_lagrangian_nodal_coordinates();
782  }
783 
784  double max_area;
785  double min_area;
786  this->max_and_min_element_size(max_area, min_area);
787  oomph_info << "Max/min element size in adapted mesh: "
788  << max_area << " "
789  << min_area << std::endl;
790  }
791  else
792  {
793  oomph_info << "Not enough benefit in adaptation.\n";
794  Nrefined=0;
795  Nunrefined=0;
796  }
797  }
798 
799 }
800 
801 #endif
refineable_tetgen_mesh.template.h
oomph::TetgenMesh::tetgenio_pt
tetgenio *& tetgenio_pt()
Access to the triangulateio representation of the mesh.
Definition: tetgen_mesh.template.h:696
oomph::RefineableTetgenMesh::adapt
void adapt(const Vector< double > &elem_error)
Adapt mesh, based on elemental error provided.
Definition: refineable_tetgen_mesh.template.cc:52