libMesh
mesh_smoother_vsmoother.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/libmesh_config.h"
19 #if defined(LIBMESH_ENABLE_VSMOOTHER)
20 
21 // Local includes
22 #include "libmesh/mesh_smoother_vsmoother.h"
23 #include "libmesh/mesh_tools.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/unstructured_mesh.h"
26 #include "libmesh/utility.h"
27 #include "libmesh/boundary_info.h"
28 #include "libmesh/equation_systems.h"
29 #include "libmesh/distributed_mesh.h"
30 #include "libmesh/steady_solver.h"
31 #include "libmesh/diff_solver.h"
32 #include "libmesh/parallel_ghost_sync.h"
33 
34 // C++ includes
35 #include <time.h> // for clock_t, clock()
36 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
37 #include <cmath>
38 #include <iomanip>
39 #include <limits>
40 
41 namespace libMesh
42 {
43 
44 // Optimization at -O2 or greater seem to break Intel's icc. So if we are
45 // being compiled with icc let's dumb-down the optimizations for this file
46 #ifdef __INTEL_COMPILER
47 # pragma optimize ( "", off )
48 #endif
49 
50 // Member functions for the Variational Smoother
52  Real dilation_weight,
53  const bool preserve_subdomain_boundaries,
54  const double relative_residual_tolerance,
55  const double absolute_residual_tolerance,
56  const bool solver_quiet,
57  const bool solver_verbose)
58  : MeshSmoother(mesh),
59  _dilation_weight(dilation_weight),
60  _preserve_subdomain_boundaries(preserve_subdomain_boundaries),
61  _setup_called(false),
62  _relative_residual_tolerance(relative_residual_tolerance),
63  _absolute_residual_tolerance(absolute_residual_tolerance),
64  _solver_quiet(solver_quiet),
65  _solver_verbose(solver_verbose)
66 {}
67 
69 {
70  // Check for multiple dimensions
71  if (_mesh.elem_dimensions().size() > 1)
72  libmesh_not_implemented_msg("Meshes containing elements of differing dimension are not yet supported.");
73 
74  /*
75  Ideally we'd want to update _mesh directly even if it already has an
76  EquationSystems attached and doing work on it ... and that's thwarted by the
77  problems:
78 
79  1. We can't easily tell if there's already an EquationSystems attached.
80  2. We can't attach a second EquationSystems safely (if it already has a system 0)
81  because the DoF indexing will need to be overwritten by our system.
82  3. The destructor of es won't even clean up after itself (we generally expect
83  a mesh to go unused after its EquationSystems is destroyed), much less know
84  how to restore anything from a previous EquationSystems.
85 
86  To avoid these issues, we'll just construct a new DistributedMesh _mesh_copy
87  from _mesh, then do the solve on _mesh_copy, then copy its node locations back
88  to _mesh after the solve is done. That'll be slightly less memory-efficient
89  and somewhat more CPU-efficient in the case where _mesh is serial, though
90  it'll be significantly less memory-efficient when _mesh is already distributed,
91  but either way the robustness is probably worth it.
92  */
93 
94  // Create a new mesh, EquationSystems, and System
95  _mesh_copy = std::make_unique<DistributedMesh>(_mesh);
96  _equation_systems = std::make_unique<EquationSystems>(*_mesh_copy);
97  _system =
98  &(_equation_systems->add_system<VariationalSmootherSystem>("variational_smoother_system"));
99 
100  // Set this to something > 0 to add more quadrature points than the default
101  // rule that integrates order 2 * fe_order + 1 polynomials exactly.
102  // Using higher quadrature orders has not had a significant effect on observed solutions.
103  //system()->extra_quadrature_order = 0;
104 
105  // Uncomment these to debug
106  //system()->print_element_solutions=true;
107  //system()->print_element_residuals=true;
108  //system()->print_element_jacobians=true;
109 
110  // Add boundary node and hanging node constraints
111  _constraint =
112  std::make_unique<VariationalSmootherConstraint>(*_system, _preserve_subdomain_boundaries);
114 
115  // Set system parameters
117 
118  // Set up solver
119  system()->time_solver = std::make_unique<SteadySolver>(*_system);
120 
121  // Uncomment this line and use -snes_test_jacobian and -snes_test_jacobian_view
122  // flags to compare the hand-coded jacobian in VariationalSmootherSystem
123  // to finite difference jacobians.
124  //system()->time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(*_system);
125 
126  _equation_systems->init();
127 
128  // Solver verbosity
129  DiffSolver & solver = *(system()->time_solver->diff_solver().get());
130  solver.quiet = _solver_quiet;
131  solver.verbose = _solver_verbose;
132 
133  // Solver convergence tolerances
134  system()->time_solver->diff_solver()->relative_residual_tolerance = _relative_residual_tolerance;
135  system()->time_solver->diff_solver()->absolute_residual_tolerance = _absolute_residual_tolerance;
136 
137  _setup_called = true;
138 }
139 
141 {
142  if (!_setup_called)
143  setup();
144 
145  system()->solve();
146 
147  // Update _mesh from _mesh_copy
148  for (auto * node_copy : _mesh_copy->local_node_ptr_range())
149  {
150  auto & node = _mesh.node_ref(node_copy->id());
151  for (const auto d : make_range(_mesh_copy->mesh_dimension()))
152  node(d) = (*node_copy)(d);
153  }
154 
155  SyncNodalPositions sync_object(_mesh);
156  Parallel::sync_dofobject_data_by_id (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
157 
158  // Release memory occupied by _mesh_copy
159  // Destruct this before _mesh_copy because it references _mesh_copy
160  _equation_systems.reset();
161  _mesh_copy.reset();
162  // We'll need to call setup again since we'll have to reconstruct _mesh_copy
163  _setup_called = false;
164 }
165 
167 {
168  libmesh_error_msg_if(!_setup_called, "Need to first call the setup() method of "
169  << "this VariationalMeshSmoother object, then call get_mesh_info() prior "
170  << "to calling smooth().");
171  return _system->get_mesh_info();
172 }
173 
174 } // namespace libMesh
175 
176 #endif // defined(LIBMESH_ENABLE_VSMOOTHER)
virtual void solve() override
Solves the system to smooth the mesh.
const Real _dilation_weight
Smoother control variables.
bool _solver_verbose
Solver verbose setting.
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:160
std::unique_ptr< EquationSystems > _equation_systems
EquationsSystems object associated with the smoother.
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:253
MeshBase & mesh
const MeshQualityInfo & get_mesh_info() const
Getter for the _system&#39;s _mesh_info attribute.
const Parallel::Communicator & comm() const
VariationalSmootherSystem * _system
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
VariationalSmootherSystem * system() const
Getter for _system to protect against dangling pointers.
The libMesh namespace provides an interface to certain functionality in the library.
std::unique_ptr< VariationalSmootherConstraint > _constraint
Constraints imposed on the smoothing process.
void attach_constraint_object(Constraint &constrain)
Register a user object for imposing constraints.
Definition: system.C:1987
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:68
std::unique_ptr< DistributedMesh > _mesh_copy
Mesh copy to avoid multiple EquationSystems.
Struct to hold smoother-relevant information about the mesh quality.
double _absolute_residual_tolerance
Solver absolute residual tolerance.
bool _setup_called
Attribute the keep track of whether the setup method has been called.
The UnstructuredMesh class is derived from the MeshBase class.
const MeshQualityInfo & get_mesh_info()
Getter for the _mesh_info attribute.
VariationalMeshSmoother(UnstructuredMesh &mesh, Real dilation_weight=0.5, const bool preserve_subdomain_boundaries=true, const double relative_residual_tolerance=TOLERANCE *TOLERANCE, const double absolute_residual_tolerance=TOLERANCE *TOLERANCE, const bool solver_quiet=true, const bool solver_verbose=false)
Constructor.
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:282
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
bool _solver_quiet
Solver quiet setting.
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:166
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
double _relative_residual_tolerance
Solver relative residual tolerance.
This class provides the necessary interface for mesh smoothing.
Definition: mesh_smoother.h:38
virtual void setup()
Setup method that creates equation systems, system, and constraints, to be called just prior to smoot...
const bool _preserve_subdomain_boundaries
Whether subdomain boundaries are subject to change via smoothing.
virtual void smooth() override
Redefinition of the smooth function from the base class.