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