Line data Source code
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 :
19 : // libmesh includes
20 : #include "libmesh/libmesh_common.h"
21 : #include "libmesh/smoothness_estimator.h"
22 : #include "libmesh/dof_map.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/error_vector.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/quadrature.h"
28 : #include "libmesh/system.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/numeric_vector.h"
31 : #include "libmesh/enum_to_string.h"
32 :
33 : // C++ includes
34 : #include <algorithm> // for std::fill
35 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
36 : #include <cmath> // for std::sqrt std::pow std::abs
37 :
38 : namespace libMesh
39 : {
40 :
41 :
42 : //-----------------------------------------------------------------
43 : // SmoothnessEstimator implementations
44 2414 : SmoothnessEstimator::SmoothnessEstimator() :
45 2414 : _extra_order(1)
46 : {
47 2414 : }
48 :
49 :
50 :
51 50172 : std::vector<Real> SmoothnessEstimator::legepoly(const unsigned int dim,
52 : const Order order,
53 : const Point p,
54 : const unsigned int matsize)
55 : {
56 4186 : std::vector<Real> psi;
57 50172 : psi.reserve(matsize);
58 :
59 : // Evaluate 1D Legendre polynomials at x, y, z up to order
60 50172 : const int npols = static_cast<int>(order) + 1;
61 :
62 54358 : std::vector<Real> Lx(npols, 0.), Ly(npols, 1.), Lz(npols, 1.);
63 50172 : const Real x = p(0);
64 50172 : Lx[0] = 1.;
65 50172 : if (npols > 1)
66 50172 : Lx[1] = x;
67 124136 : for (int n = 2; n < npols; ++n)
68 92504 : Lx[n] = ((2. * n - 1) * x * Lx[n - 1] - (n - 1) * Lx[n - 2]) / n;
69 :
70 50172 : if (dim > 1)
71 : {
72 40320 : const Real y = p(1);
73 40320 : Ly[0] = 1.;
74 40320 : if (npols > 1)
75 40320 : Ly[1] = y;
76 80640 : for (int n = 2; n < npols; ++n)
77 50400 : Ly[n] = ((2. * n - 1) * y * Ly[n - 1] - (n - 1) * Ly[n - 2]) / n;
78 : }
79 :
80 41972 : if (dim > 2)
81 : {
82 0 : const Real z = p(2);
83 0 : Lz[0] = 1.;
84 0 : if (npols > 1)
85 0 : Lz[1] = z;
86 0 : for (int n = 2; n < npols; ++n)
87 0 : Lz[n] = ((2. * n - 1) * z * Lz[n - 1] - (n - 1) * Lz[n - 2]) / n;
88 : }
89 :
90 : // Loop over total degree and build tensor-product Legendre basis
91 224480 : for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(order); poly_deg++)
92 : {
93 174308 : switch (dim)
94 : {
95 0 : case 3:
96 0 : for (int i = poly_deg; i >= 0; --i)
97 0 : for (int j = poly_deg - i; j >= 0; --j)
98 : {
99 0 : int k = poly_deg - i - j;
100 0 : psi.push_back(Lx[i] * Ly[j] * Lz[k]);
101 0 : }
102 0 : break;
103 :
104 120960 : case 2:
105 362880 : for (int i = poly_deg; i >= 0; --i)
106 : {
107 241920 : int j = poly_deg - i;
108 282240 : psi.push_back(Lx[i] * Ly[j]);
109 20160 : }
110 10080 : break;
111 :
112 53348 : case 1:
113 57820 : psi.push_back(Lx[poly_deg]);
114 4472 : break;
115 :
116 0 : default:
117 0 : libmesh_error_msg("Invalid dimension dim " << dim);
118 : }
119 : }
120 :
121 54358 : return psi;
122 : }
123 :
124 4166 : Real SmoothnessEstimator::compute_slope(int N, Real Sx, Real Sy, Real Sxx, Real Sxy)
125 : {
126 4166 : const Real denom = (N * Sxx - Sx * Sx);
127 : // Triggers for first order polynomials when there only 1 point
128 : // available at log |c_k| and log |k| space to fit.
129 4166 : if (std::abs(denom) < std::numeric_limits<Real>::epsilon()) {
130 0 : return std::numeric_limits<Real>::max();
131 : }
132 4166 : return (N * Sxy - Sx * Sy) / denom;
133 : }
134 :
135 2414 : void SmoothnessEstimator::reduce_smoothness (std::vector<ErrorVectorReal> & smoothness_per_cell,
136 : const Parallel::Communicator & comm)
137 : {
138 : // Aggregates element-wise contributions computed
139 : // in parallel across all processors
140 2414 : comm.sum(smoothness_per_cell);
141 2414 : }
142 :
143 2414 : void SmoothnessEstimator::estimate_smoothness (const System & system,
144 : ErrorVector & smoothness_per_cell,
145 : const NumericVector<Number> * solution_vector)
146 : {
147 136 : LOG_SCOPE("estimate_smoothness()", "SmoothnessEstimator");
148 :
149 : // The current mesh
150 136 : const MeshBase & mesh = system.get_mesh();
151 :
152 : // Resize the smoothness_per_cell vector to be
153 : // the number of elements, initialize it to 0.
154 2414 : smoothness_per_cell.resize (mesh.max_elem_id());
155 68 : std::fill (smoothness_per_cell.begin(), smoothness_per_cell.end(), 0.);
156 :
157 : // Prepare current_local_solution to localize a non-standard
158 : // solution vector if necessary
159 2414 : if (solution_vector && solution_vector != system.solution.get())
160 : {
161 0 : NumericVector<Number> * newsol =
162 : const_cast<NumericVector<Number> *>(solution_vector);
163 0 : System & sys = const_cast<System &>(system);
164 0 : newsol->swap(*sys.solution);
165 0 : sys.update();
166 : }
167 :
168 : //------------------------------------------------------------
169 : // Iterate over all the active elements in the mesh
170 : // that live on this processor.
171 4760 : Threads::parallel_for (ConstElemRange(mesh.active_local_elements_begin(),
172 2482 : mesh.active_local_elements_end(),
173 : 200),
174 2346 : EstimateSmoothness(system,
175 : *this,
176 : smoothness_per_cell)
177 : );
178 :
179 : // Each processor has now computed the smoothness for its local
180 : // elements, and smoothness_per_cell contains 0 for all the
181 : // non-local elements. Summing the vector will provide the true
182 : // value for each element, local or remote
183 2414 : this->reduce_smoothness (smoothness_per_cell, system.comm());
184 :
185 : // If we used a non-standard solution before, now is the time to fix
186 : // the current_local_solution
187 2414 : if (solution_vector && solution_vector != system.solution.get())
188 : {
189 0 : NumericVector<Number> * newsol = const_cast<NumericVector<Number> *>(solution_vector);
190 0 : System & sys = const_cast<System &>(system);
191 0 : newsol->swap(*sys.solution);
192 0 : sys.update();
193 : }
194 2414 : }
195 :
196 :
197 2552 : void SmoothnessEstimator::EstimateSmoothness::operator()(const ConstElemRange & range) const
198 : {
199 : // The current mesh
200 2552 : const MeshBase & mesh = system.get_mesh();
201 :
202 : // The dimensionality of the mesh
203 2552 : const unsigned int dim = mesh.mesh_dimension();
204 :
205 : // The number of variables in the system
206 2552 : const unsigned int n_vars = system.n_vars();
207 :
208 : // The DofMap for this system
209 114 : const DofMap & dof_map = system.get_dof_map();
210 :
211 : //------------------------------------------------------------
212 : // Iterate over all the elements in the range.
213 6718 : for (const auto & elem : range)
214 : {
215 4166 : const dof_id_type e_id = elem->id();
216 : // Loop over each variable
217 8332 : for (unsigned int var=0; var<n_vars; var++)
218 : {
219 348 : const FEType & fe_type = dof_map.variable_type(var);
220 4166 : const Order element_order = fe_type.order + elem->p_level();
221 :
222 4514 : std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
223 4514 : std::unique_ptr<QBase> qrule(fe_type.default_quadrature_rule(dim, smoothness_estimator._extra_order));
224 4514 : fe->attach_quadrature_rule(qrule.get());
225 :
226 1044 : const std::vector<Real> & JxW = fe->get_JxW();
227 1044 : const std::vector<Point> & q_point = fe->get_xyz();
228 348 : const std::vector<std::vector<Real>> * phi = &(fe->get_phi());
229 :
230 696 : std::vector<dof_id_type> dof_indices;
231 :
232 4166 : unsigned int matsize = element_order + 1;
233 4166 : if (dim > 1)
234 : {
235 2520 : matsize *= (element_order + 2);
236 2520 : matsize /= 2;
237 : }
238 2796 : if (dim > 2)
239 : {
240 0 : matsize *= (element_order + 3);
241 0 : matsize /= 3;
242 : }
243 :
244 4862 : DenseMatrix<Number> Kp(matsize, matsize);
245 4166 : DenseVector<Number> F;
246 4166 : DenseVector<Number> Pu_h;
247 :
248 3818 : F.resize(matsize);
249 3818 : Pu_h.resize(matsize);
250 :
251 :
252 4166 : fe->reinit(elem);
253 :
254 4166 : dof_map.dof_indices(elem, dof_indices, var);
255 348 : libmesh_assert_equal_to(dof_indices.size(), phi->size());
256 :
257 696 : const unsigned int n_dofs = cast_int<unsigned int>(dof_indices.size());
258 348 : const unsigned int n_qp = qrule->n_points();
259 :
260 54338 : for (unsigned int qp = 0; qp < n_qp; qp++)
261 : {
262 58544 : std::vector<Real> psi = legepoly(dim, element_order, q_point[qp], matsize);
263 8372 : const unsigned int psi_size = cast_int<unsigned int>(psi.size());
264 :
265 345440 : for (unsigned int i = 0; i < matsize; i++)
266 2065504 : for (unsigned int j = 0; j < matsize; j++)
267 2360908 : Kp(i, j) += JxW[qp] * psi[i] * psi[j];
268 :
269 4186 : Number u_h = libMesh::zero;
270 466400 : for (unsigned int i = 0; i < n_dofs; i++)
271 520364 : u_h += (*phi)[i][qp] * system.current_solution(dof_indices[i]);
272 :
273 345440 : for (unsigned int i = 0; i < psi_size; i++)
274 369164 : F(i) += JxW[qp] * u_h * psi[i];
275 : }
276 :
277 4166 : Kp.lu_solve(F, Pu_h);
278 :
279 : // Generate index to total degree map. Total_degree_per_index[i] gives the degree ith basis function
280 696 : std::vector<unsigned int> total_degree_per_index;
281 :
282 19932 : for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(element_order); ++poly_deg)
283 : {
284 15766 : switch (dim)
285 : {
286 0 : case 3:
287 0 : for (int i = poly_deg; i >= 0; --i)
288 0 : for (int j = poly_deg - i; j >= 0; --j)
289 : {
290 0 : int k = poly_deg - i - j;
291 0 : total_degree_per_index.push_back(i + j + k);
292 0 : }
293 0 : break;
294 :
295 7560 : case 2:
296 22680 : for (int i = poly_deg; i >= 0; --i)
297 : {
298 2520 : int j = poly_deg - i;
299 15120 : total_degree_per_index.push_back(i + j);
300 1260 : }
301 630 : break;
302 :
303 8206 : case 1:
304 8206 : total_degree_per_index.push_back(poly_deg);
305 688 : break;
306 :
307 0 : default:
308 0 : libmesh_error_msg("Invalid dimension dim " << dim);
309 : }
310 : }
311 :
312 : // Group coefficients |c_k| by degree k
313 696 : std::map<unsigned int, std::vector<Number>> coeff_by_degree;
314 :
315 27492 : for (unsigned int i = 0; i < Pu_h.size(); ++i)
316 : {
317 23326 : unsigned int degree = total_degree_per_index[i];
318 : // Excluding the constant mode i.e, zeroth order coefficient as we use ln|order| later
319 : // Constant order term often dominates the scale and doesn't inform about regularity.
320 23326 : if (degree > 0)
321 19160 : coeff_by_degree[degree].push_back(Pu_h(i));
322 : }
323 :
324 : // Compute L2 norm of each group |c_k|
325 696 : std::vector<unsigned int> degrees;
326 696 : std::vector<double> norms;
327 :
328 15766 : for (const auto & pair : coeff_by_degree)
329 : {
330 11600 : unsigned int deg = pair.first;
331 970 : const std::vector<Number> & coeffs = pair.second;
332 :
333 11600 : double norm = 0.;
334 30760 : for (const auto & c : coeffs)
335 20760 : norm += std::norm(c);
336 :
337 11600 : norm = std::sqrt(norm);
338 :
339 11600 : degrees.push_back(deg);
340 11600 : norms.push_back(norm);
341 : }
342 :
343 :
344 : // Collect log |c_k| and log |k|
345 696 : std::vector<double> log_norms, log_degrees;
346 :
347 16114 : for (unsigned int i = 0; i < degrees.size(); ++i)
348 : {
349 : // filter tiny terms
350 12570 : if (norms[i] > 1e-12)
351 : {
352 11600 : log_degrees.push_back(std::log(static_cast<double>(degrees[i])));
353 12570 : log_norms.push_back(std::log(norms[i]));
354 : }
355 : }
356 :
357 :
358 : // Fit log-log line - we use least-squares fit
359 348 : double Sx = 0., Sy = 0., Sxx = 0., Sxy = 0.;
360 696 : const size_t N = log_degrees.size();
361 :
362 15766 : for (size_t i = 0; i < N; ++i)
363 : {
364 11600 : Sx += log_degrees[i];
365 11600 : Sy += log_norms[i];
366 11600 : Sxx += log_degrees[i] * log_degrees[i];
367 11600 : Sxy += log_degrees[i] * log_norms[i];
368 : }
369 :
370 4166 : const double regularity = -compute_slope(N, Sx, Sy, Sxx, Sxy);
371 : // const double intercept = (Sy - slope * Sx) / N;
372 :
373 348 : Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
374 4514 : smoothness_per_cell[e_id] = regularity;
375 6940 : } // end variables loop
376 :
377 : } // end element loop
378 :
379 2552 : } // End () operator definition
380 :
381 : } // namespace libMesh
|
