gmsh_tet_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_GMSH_TET_MESH_TEMPLATE_CC
32 #define OOMPH_GMSH_TET_MESH_TEMPLATE_CC
33 
34 
35 #include "gmsh_tet_mesh.template.h"
36 
37 
38 namespace oomph
39 {
40 
41 //======================================================================
42 /// Adapt problem based on specified elemental error estimates
43 //======================================================================
44 template <class ELEMENT>
45 void RefineableGmshTetMesh<ELEMENT>::adapt(const Vector<double>& elem_error)
46  {
47 
48 
49  double t_start=0.0;
50 
51  //###################################
52  t_start=TimingHelpers::timer();
53  //###################################
54 
55  // Get refinement targets
56  Vector<double> target_size(elem_error.size());
57  double max_edge_ratio=this->compute_volume_target(elem_error,
58  target_size);
59  // Get maximum target volume
60  unsigned n=target_size.size();
61  double max_size=0.0;
62  double min_size=DBL_MAX;
63  for (unsigned e=0;e<n;e++)
64  {
65  if (target_size[e]>max_size) max_size=target_size[e];
66  if (target_size[e]<min_size) min_size=target_size[e];
67  }
68 
69  oomph_info << "Maximum target size: " << max_size << std::endl;
70  oomph_info << "Minimum target size: " << min_size << std::endl;
71  oomph_info << "Number of elements to be refined "
72  << this->Nrefined << std::endl;
73  oomph_info << "Number of elements to be unrefined "
74  << this->Nunrefined << std::endl;
75  oomph_info << "Max edge ratio "<< max_edge_ratio << std::endl;
76 
77  double orig_max_size, orig_min_size;
78  this->max_and_min_element_size(orig_max_size, orig_min_size);
79  oomph_info << "Max/min element size in original mesh: "
80  << orig_max_size << " "
81  << orig_min_size << std::endl;
82 
83  //##################################################################
84  oomph_info << "adapt: Time for getting volume targets : "
85  << TimingHelpers::timer()-t_start
86  << " sec " << std::endl;
87  //##################################################################
88 
89  // Should we bother to adapt?
90  if ( (Nrefined > 0) || (Nunrefined > this->max_keep_unrefined())
91  || (max_edge_ratio > this->max_permitted_edge_ratio()) )
92  {
93 
94  if (! ( (Nrefined > 0) || (Nunrefined > max_keep_unrefined()) ) )
95  {
96  oomph_info
97  << "Mesh regeneration triggered by edge ratio criterion\n";
98  }
99 
100 
101  // Are we dealing with a solid mesh?
102  SolidMesh* solid_mesh_pt=dynamic_cast<SolidMesh*>(this);
103 
104 
105  // If the mesh is a solid mesh then do the mapping based on the
106  // Eulerian coordinates
107  bool use_eulerian_coords=false;
108  if (solid_mesh_pt!=0)
109  {
110  use_eulerian_coords=true;
111  }
112 
113  MeshAsGeomObject* mesh_geom_obj_pt=0;
114 
115 #ifdef OOMPH_HAS_CGAL
116 
117  // Make cgal-based bin
118  CGALSamplePointContainerParameters cgal_params(this);
119  if (use_eulerian_coords)
120  {
121  cgal_params.enable_use_eulerian_coordinates_during_setup();
122  }
123  mesh_geom_obj_pt=new MeshAsGeomObject(&cgal_params);
124 
125 #else
126 
127  std::ostringstream error_message;
128  error_message << "Non-CGAL-based target size transfer not implemented.\n";
129  throw OomphLibError(error_message.str(),
130  OOMPH_CURRENT_FUNCTION,
131  OOMPH_EXCEPTION_LOCATION);
132 
133  // Do something here...
134 
135 #endif
136 
137  // Set up a map from pointer to element to its number
138  // in the mesh
139  std::map<GeneralisedElement*,unsigned> element_number;
140  unsigned nelem=this->nelement();
141  for (unsigned e=0;e<nelem;e++)
142  {
143  element_number[this->element_pt(e)]=e;
144  }
145 
146  // Get min/max coordinates
147  Vector<std::pair<double, double> > min_and_max_coordinates=
148  mesh_geom_obj_pt->sample_point_container_pt()->
149  min_and_max_coordinates();
150 
151  // Setup grid dimensions for size transfer
152  double dx=(min_and_max_coordinates[0].second-
153  min_and_max_coordinates[0].first);
154  double dy=(min_and_max_coordinates[1].second-
155  min_and_max_coordinates[1].first);
156  double dz=(min_and_max_coordinates[2].second-
157  min_and_max_coordinates[2].first);
158 
159  double dx_target=this->Gmsh_parameters_pt->dx_for_target_size_transfer();
160  unsigned nx=unsigned(dx/dx_target);
161  unsigned ny=unsigned(dy/dx_target);
162  unsigned nz=unsigned(dz/dx_target);
163 
164  dx/=double(nx);
165  dy/=double(ny);
166  dz/=double(nz);
167 
168 
169  // Size transfer via hard disk -- yikes...
170  std::string target_size_file_name=
171  this->Gmsh_parameters_pt->target_size_file_name();
172 
173  std::ofstream target_size_file;
174  target_size_file.open(target_size_file_name.c_str());
175  target_size_file << min_and_max_coordinates[0].first << " "
176  << min_and_max_coordinates[1].first << " "
177  << min_and_max_coordinates[2].first << " "
178  << std::endl;
179  target_size_file << dx << " " << dy << " " << dz << std::endl;
180  target_size_file << nx+1 << " " << ny+1 << " " << nz+1 << std::endl;
181 
182 
183  // Doc target areas
184  int counter=this->Gmsh_parameters_pt->
185  counter_for_filename_gmsh_size_transfer();
186  std::string stem=this->Gmsh_parameters_pt->
187  stem_for_filename_gmsh_size_transfer();
188  std::ofstream latest_sizes_file;
189  bool doc_target_areas=false;
190  if ((stem!="")&&(counter>=0))
191  {
192  doc_target_areas=true;
193  std::string filename=stem+
194  oomph::StringConversion::to_string(counter)+".dat";
195  latest_sizes_file.open(filename.c_str());
196  latest_sizes_file << "ZONE I=" << nx+1
197  << ", J=" << ny+1
198  << ", K=" << nz+1
199  << std::endl;
200  this->Gmsh_parameters_pt->
201  counter_for_filename_gmsh_size_transfer()++;
202  }
203 
204 
205  Vector<double> x(3);
206  for (unsigned i=0;i<=nx;i++)
207  {
208  x[0]=min_and_max_coordinates[0].first+double(i)*dx;
209  for (unsigned j=0;j<=ny;j++)
210  {
211  x[1]=min_and_max_coordinates[1].first+double(j)*dy;
212  for (unsigned k=0;k<=nz;k++)
213  {
214  x[2]=min_and_max_coordinates[2].first+double(k)*dz;
215 
216  // Try the specified number of nearest sample points for Newton search
217  // then just settle on the nearest one
218  Vector<double> s(3);
219  GeomObject* geom_obj_pt=0;
220  unsigned max_sample_points=this->Gmsh_parameters_pt->
221  max_sample_points_for_limited_locate_zeta_during_target_size_transfer();
222 
223 #ifdef OOMPH_HAS_CGAL
224 
225  dynamic_cast<CGALSamplePointContainer*>(mesh_geom_obj_pt->
226  sample_point_container_pt())->
227  limited_locate_zeta(x,max_sample_points,
228  geom_obj_pt,s);
229 
230 #else
231 
232  std::ostringstream error_message;
233  error_message << "Non-CGAL-based target size transfer not implemented.\n";
234  throw OomphLibError(error_message.str(),
235  OOMPH_CURRENT_FUNCTION,
236  OOMPH_EXCEPTION_LOCATION);
237 
238  // Do something here...
239 
240 #endif
241 
242 
243 #ifdef PARANOID
244  if (geom_obj_pt==0)
245  {
246  std::stringstream error_message;
247  error_message
248  << "Limited locate zeta failed for zeta = [ "
249  << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
250  throw OomphLibError(error_message.str(),
251  OOMPH_CURRENT_FUNCTION,
252  OOMPH_EXCEPTION_LOCATION);
253  }
254  else
255  {
256 #endif
257 
258  FiniteElement* fe_pt=dynamic_cast<FiniteElement*>(geom_obj_pt);
259 
260 #ifdef PARANOID
261  if (fe_pt==0)
262  {
263  std::stringstream error_message;
264  error_message
265  << "Cast to FE for GeomObject returned by limited "
266  << "locate zeta failed for zeta = [ "
267  << x[0] << " " << x[1] << " " << x[2] << " ]. Makes no sense!\n";
268  throw OomphLibError(error_message.str(),
269  OOMPH_CURRENT_FUNCTION,
270  OOMPH_EXCEPTION_LOCATION);
271  }
272  else
273  {
274 #endif
275 
276  // What's the target size of the element that contains this point
277  double tg_size=pow(target_size[element_number[fe_pt]],1.0/3.0);
278  target_size_file << tg_size << " ";
279 
280  // Doc?
281  if (doc_target_areas)
282  {
283  latest_sizes_file << x[0] << " "
284  << x[1] << " "
285  << x[2] << " "
286  << tg_size << std::endl;
287  }
288 
289 #ifdef PARANOID
290  }
291  }
292 #endif
293  }
294  target_size_file << std::endl;
295  }
296  }
297  target_size_file.close();
298 
299  if (doc_target_areas)
300  {
301  latest_sizes_file.close();
302  }
303 
304  // Build new mesh, reading area information from file
305  bool use_mesh_grading_from_file=true;
306  RefineableGmshTetMesh<ELEMENT>* new_mesh_pt= new
307  RefineableGmshTetMesh<ELEMENT>(this->Gmsh_parameters_pt,
308  use_mesh_grading_from_file,
309  this->Time_stepper_pt);
310 
311  /* // keep around for debugging */
312  /* unsigned nplot=5; */
313  /* ofstream vtu_file; */
314  /* vtu_file.open("new_mesh.vtu"); */
315  /* new_mesh_pt->output_paraview(vtu_file,nplot); */
316  /* vtu_file.close(); */
317 
318  //###################################
319  t_start=TimingHelpers::timer();
320  //###################################
321 
322  ProjectionProblem<ELEMENT>* project_problem_pt=0;
323 
324  // Check that the projection step is not disabled
325  if (!this->Gmsh_parameters_pt->projection_is_disabled())
326  {
327  // Project current solution onto new mesh
328  //---------------------------------------
329  project_problem_pt=new ProjectionProblem<ELEMENT>;
330  project_problem_pt->mesh_pt()=new_mesh_pt;
331  project_problem_pt->project(this);
332 
333  oomph_info
334  << "adapt: Time for project soln onto new mesh : "
335  << TimingHelpers::timer()-t_start
336  << " sec " << std::endl;
337  }
338  //###################################
339  t_start=TimingHelpers::timer();
340  //###################################
341 
342  //Flush the old mesh
343  unsigned nnod=nnode();
344  for(unsigned j=nnod;j>0;j--)
345  {
346  delete Node_pt[j-1];
347  Node_pt[j-1] = 0;
348  }
349  unsigned nel=nelement();
350  for(unsigned e=nel;e>0;e--)
351  {
352  delete Element_pt[e-1];
353  Element_pt[e-1] = 0;
354  }
355 
356  // Now copy back to current mesh
357  //------------------------------
358  nnod=new_mesh_pt->nnode();
359  Node_pt.resize(nnod);
360  nel=new_mesh_pt->nelement();
361  Element_pt.resize(nel);
362  for(unsigned j=0;j<nnod;j++)
363  {
364  Node_pt[j] = new_mesh_pt->node_pt(j);
365  }
366  for(unsigned e=0;e<nel;e++)
367  {
368  Element_pt[e] = new_mesh_pt->element_pt(e);
369  }
370 
371  //Copy the boundary schemes
372  unsigned nbound=new_mesh_pt->nboundary();
373  Boundary_element_pt.resize(nbound);
374  Face_index_at_boundary.resize(nbound);
375  Boundary_node_pt.resize(nbound);
376  for (unsigned b=0;b<nbound;b++)
377  {
378  unsigned nel=new_mesh_pt->nboundary_element(b);
379  Boundary_element_pt[b].resize(nel);
380  Face_index_at_boundary[b].resize(nel);
381  for (unsigned e=0;e<nel;e++)
382  {
383  Boundary_element_pt[b][e]=new_mesh_pt->boundary_element_pt(b,e);
384  Face_index_at_boundary[b][e]=new_mesh_pt->face_index_at_boundary(b,e);
385  }
386  unsigned nnod=new_mesh_pt->nboundary_node(b);
387  Boundary_node_pt[b].resize(nnod);
388  for (unsigned j=0;j<nnod;j++)
389  {
390  Boundary_node_pt[b][j]=new_mesh_pt->boundary_node_pt(b,j);
391  }
392  }
393 
394  //Also copy over the new boundary and region information
395  unsigned n_region = new_mesh_pt->nregion();
396 
397  //Only bother if we have regions
398  if(n_region > 1)
399  {
400  //Deal with the region information first
401  this->Region_element_pt.resize(n_region);
402  this->Region_attribute.resize(n_region);
403  for(unsigned i=0;i<n_region;i++)
404  {
405  // Copy across region attributes (region ids!)
406  this->Region_attribute[i] = new_mesh_pt->region_attribute(i);
407 
408  //Find the number of elements in the region
409  unsigned r=this->Region_attribute[i];
410  unsigned n_region_element = new_mesh_pt->nregion_element(r);
411  this->Region_element_pt[i].resize(n_region_element);
412  for(unsigned e=0;e<n_region_element;e++)
413  {
414  this->Region_element_pt[i][e] = new_mesh_pt->region_element_pt(r,e);
415  }
416  }
417 
418  //Now the boundary region information
419  this->Boundary_region_element_pt.resize(nbound);
420  this->Face_index_region_at_boundary.resize(nbound);
421 
422  //Now loop over the boundaries
423  for(unsigned b=0;b<nbound;++b)
424  {
425  //Loop over the regions
426  for(unsigned i=0;i<n_region;++i)
427  {
428  unsigned r=this->Region_attribute[i];
429 
430  unsigned n_boundary_el_in_region =
431  new_mesh_pt->nboundary_element_in_region(b,r);
432  if(n_boundary_el_in_region > 0)
433  {
434  //Allocate storage in the map
435  this->Boundary_region_element_pt[b][r].
436  resize(n_boundary_el_in_region);
437  this->Face_index_region_at_boundary[b][r].
438  resize(n_boundary_el_in_region);
439 
440  //Copy over the information
441  for(unsigned e=0;e<n_boundary_el_in_region;++e)
442  {
443  this->Boundary_region_element_pt[b][r][e]
444  = new_mesh_pt->boundary_element_in_region_pt(b,r,e);
445  this->Face_index_region_at_boundary[b][r][e]
446  = new_mesh_pt->face_index_at_boundary_in_region(b,r,e);
447  }
448  }
449  }
450  } //End of loop over boundaries
451 
452  } //End of case when more than one region
453 
454 
455  // Flush the mesh
456  new_mesh_pt->flush_element_and_node_storage();
457 
458  // Delete the mesh and the problem
459  delete new_mesh_pt;
460  delete project_problem_pt;
461 
462  //##################################################################
463  oomph_info << "adapt: Time for moving nodes etc. to actual mesh : "
464  << TimingHelpers::timer()-t_start
465  << " sec " << std::endl;
466  //##################################################################
467 
468 
469  // Solid mesh?
470  if (solid_mesh_pt!=0)
471  {
472  // Warning
473  std::stringstream error_message;
474  error_message
475  << "Lagrangian coordinates are currently not projected but are\n"
476  << "are re-set during adaptation. This is not appropriate for\n"
477  << "real solid mechanics problems!\n";
478  OomphLibWarning(error_message.str(),
479  OOMPH_CURRENT_FUNCTION,
480  OOMPH_EXCEPTION_LOCATION);
481 
482  // Reset Lagrangian coordinates
483  dynamic_cast<SolidMesh*>(this)->set_lagrangian_nodal_coordinates();
484  }
485 
486  double max_area;
487  double min_area;
488  this->max_and_min_element_size(max_area, min_area);
489  oomph_info << "Max/min element size in adapted mesh: "
490  << max_area << " "
491  << min_area << std::endl;
492  }
493  else
494  {
495  oomph_info << "Not enough benefit in adaptation.\n";
496  Nrefined=0;
497  Nunrefined=0;
498  }
499 }
500 }
501 
502 
503 #endif
gmsh_tet_mesh.template.h
oomph::RefineableGmshTetMesh::adapt
void adapt(const Vector< double > &elem_error)
Adapt mesh, based on elemental error provided.
Definition: gmsh_tet_mesh.template.cc:45