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  : MeshSmoother(mesh),
55  _dilation_weight(dilation_weight),
56  _preserve_subdomain_boundaries(preserve_subdomain_boundaries),
57  _setup_called(false)
58 {}
59 
61 {
62  // Check for multiple dimensions
63  if (_mesh.elem_dimensions().size() > 1)
64  libmesh_not_implemented_msg("Meshes containing elements of differing dimension are not yet supported.");
65 
66  /*
67  Ideally we'd want to update _mesh directly even if it already has an
68  EquationSystems attached and doing work on it ... and that's thwarted by the
69  problems:
70 
71  1. We can't easily tell if there's already an EquationSystems attached.
72  2. We can't attach a second EquationSystems safely (if it already has a system 0)
73  because the DoF indexing will need to be overwritten by our system.
74  3. The destructor of es won't even clean up after itself (we generally expect
75  a mesh to go unused after its EquationSystems is destroyed), much less know
76  how to restore anything from a previous EquationSystems.
77 
78  To avoid these issues, we'll just construct a new DistributedMesh _mesh_copy
79  from _mesh, then do the solve on _mesh_copy, then copy its node locations back
80  to _mesh after the solve is done. That'll be slightly less memory-efficient
81  and somewhat more CPU-efficient in the case where _mesh is serial, though
82  it'll be significantly less memory-efficient when _mesh is already distributed,
83  but either way the robustness is probably worth it.
84  */
85 
86  // Create a new mesh, EquationSystems, and System
87  _mesh_copy = std::make_unique<DistributedMesh>(_mesh);
88  _equation_systems = std::make_unique<EquationSystems>(*_mesh_copy);
89  _system =
90  &(_equation_systems->add_system<VariationalSmootherSystem>("variational_smoother_system"));
91 
92  // Set this to something > 0 to add more quadrature points than the default
93  // rule that integrates order 2 * fe_order + 1 polynomials exactly.
94  // Using higher quadrature orders has not had a significant effect on observed solutions.
95  //system()->extra_quadrature_order = 0;
96 
97  // Uncomment these to debug
98  //system()->print_element_solutions=true;
99  //system()->print_element_residuals=true;
100  //system()->print_element_jacobians=true;
101 
102  // Add boundary node and hanging node constraints
103  _constraint =
104  std::make_unique<VariationalSmootherConstraint>(*_system, _preserve_subdomain_boundaries);
106 
107  // Set system parameters
109 
110  // Set up solver
111  system()->time_solver = std::make_unique<SteadySolver>(*_system);
112 
113  // Uncomment this line and use -snes_test_jacobian and -snes_test_jacobian_view
114  // flags to compare the hand-coded jacobian in VariationalSmootherSystem
115  // to finite difference jacobians.
116  //system()->time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(*_system);
117 
118  _equation_systems->init();
119 
120  // More debugging options
121  // DiffSolver & solver = *(system()->time_solver->diff_solver().get());
122  // solver.quiet = false;
123  // solver.verbose = true;
124 
125  system()->time_solver->diff_solver()->relative_residual_tolerance = TOLERANCE * TOLERANCE;
126  system()->time_solver->diff_solver()->absolute_residual_tolerance = TOLERANCE * TOLERANCE;
127 
128  _setup_called = true;
129 }
130 
132 {
133  if (!_setup_called)
134  setup();
135 
136  system()->solve();
137 
138  // Update _mesh from _mesh_copy
139  for (auto * node_copy : _mesh_copy->local_node_ptr_range())
140  {
141  auto & node = _mesh.node_ref(node_copy->id());
142  for (const auto d : make_range(_mesh_copy->mesh_dimension()))
143  node(d) = (*node_copy)(d);
144  }
145 
146  SyncNodalPositions sync_object(_mesh);
147  Parallel::sync_dofobject_data_by_id (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
148 
149  // Release memory occupied by _mesh_copy
150  // Destruct this before _mesh_copy because it references _mesh_copy
151  _equation_systems.reset();
152  _mesh_copy.reset();
153  // We'll need to call setup again since we'll have to reconstruct _mesh_copy
154  _setup_called = false;
155 }
156 
158 {
159  libmesh_error_msg_if(!_setup_called, "Need to first call the setup() method of "
160  << "this VariationalMeshSmoother object, then call get_mesh_info() prior "
161  << "to calling smooth().");
162  return _system->get_mesh_info();
163 }
164 
165 } // namespace libMesh
166 
167 #endif // defined(LIBMESH_ENABLE_VSMOOTHER)
virtual void solve() override
Solves the system to smooth the mesh.
const Real _dilation_weight
Smoother control variables.
std::unique_ptr< EquationSystems > _equation_systems
EquationsSystems object associated with the smoother.
static constexpr Real TOLERANCE
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:2209
std::unique_ptr< DistributedMesh > _mesh_copy
Mesh copy to avoid multiple EquationSystems.
Struct to hold smoother-relevant information about the mesh quality.
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.
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.
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
VariationalMeshSmoother(UnstructuredMesh &mesh, Real dilation_weight=0.5, const bool preserve_subdomain_boundaries=true)
Simple constructor to use for smoothing purposes.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
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.