Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 50176 : 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 50176 : psi.reserve(matsize);
58 :
59 : // Evaluate 1D Legendre polynomials at x, y, z up to order
60 50176 : const int npols = static_cast<int>(order) + 1;
61 :
62 54362 : std::vector<Real> Lx(npols, 0.), Ly(npols, 1.), Lz(npols, 1.);
63 50176 : const Real x = p(0);
64 50176 : Lx[0] = 1.;
65 50176 : if (npols > 1)
66 50176 : Lx[1] = x;
67 124180 : for (int n = 2; n < npols; ++n)
68 92544 : Lx[n] = ((2. * n - 1) * x * Lx[n - 1] - (n - 1) * Lx[n - 2]) / n;
69 :
70 50176 : 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 224532 : for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(order); poly_deg++)
92 : {
93 174356 : 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 53396 : case 1:
113 57868 : 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 54362 : 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 204 : 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 2482 : Threads::parallel_for (mesh.active_local_element_stored_range(),
172 2346 : EstimateSmoothness(system,
173 : *this,
174 : smoothness_per_cell)
175 : );
176 :
177 : // Each processor has now computed the smoothness for its local
178 : // elements, and smoothness_per_cell contains 0 for all the
179 : // non-local elements. Summing the vector will provide the true
180 : // value for each element, local or remote
181 2414 : this->reduce_smoothness (smoothness_per_cell, system.comm());
182 :
183 : // If we used a non-standard solution before, now is the time to fix
184 : // the current_local_solution
185 2414 : if (solution_vector && solution_vector != system.solution.get())
186 : {
187 0 : NumericVector<Number> * newsol = const_cast<NumericVector<Number> *>(solution_vector);
188 0 : System & sys = const_cast<System &>(system);
189 0 : newsol->swap(*sys.solution);
190 0 : sys.update();
191 : }
192 2414 : }
193 :
194 :
195 2414 : void SmoothnessEstimator::EstimateSmoothness::operator()(const ConstElemRange & range) const
196 : {
197 : // The current mesh
198 2414 : const MeshBase & mesh = system.get_mesh();
199 :
200 : // The dimensionality of the mesh
201 2414 : const unsigned int dim = mesh.mesh_dimension();
202 :
203 : // The number of variables in the system
204 2414 : const unsigned int n_vars = system.n_vars();
205 :
206 : // The DofMap for this system
207 2414 : const DofMap & dof_map = system.get_dof_map();
208 :
209 : //------------------------------------------------------------
210 : // Iterate over all the elements in the range.
211 6580 : for (const auto & elem : range)
212 : {
213 4166 : const dof_id_type e_id = elem->id();
214 : // Loop over each variable
215 8332 : for (unsigned int var=0; var<n_vars; var++)
216 : {
217 348 : const FEType & fe_type = dof_map.variable_type(var);
218 4166 : const Order element_order = fe_type.order + elem->p_level();
219 :
220 4514 : std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
221 4514 : std::unique_ptr<QBase> qrule(fe_type.default_quadrature_rule(dim, smoothness_estimator._extra_order));
222 4514 : fe->attach_quadrature_rule(qrule.get());
223 :
224 1044 : const std::vector<Real> & JxW = fe->get_JxW();
225 1044 : const std::vector<Point> & q_point = fe->get_xyz();
226 348 : const std::vector<std::vector<Real>> * phi = &(fe->get_phi());
227 :
228 696 : std::vector<dof_id_type> dof_indices;
229 :
230 4166 : unsigned int matsize = element_order + 1;
231 4166 : if (dim > 1)
232 : {
233 2520 : matsize *= (element_order + 2);
234 2520 : matsize /= 2;
235 : }
236 2796 : if (dim > 2)
237 : {
238 0 : matsize *= (element_order + 3);
239 0 : matsize /= 3;
240 : }
241 :
242 4862 : DenseMatrix<Number> Kp(matsize, matsize);
243 4166 : DenseVector<Number> F;
244 4166 : DenseVector<Number> Pu_h;
245 :
246 3818 : F.resize(matsize);
247 3818 : Pu_h.resize(matsize);
248 :
249 :
250 4166 : fe->reinit(elem);
251 :
252 4166 : dof_map.dof_indices(elem, dof_indices, var);
253 348 : libmesh_assert_equal_to(dof_indices.size(), phi->size());
254 :
255 696 : const unsigned int n_dofs = cast_int<unsigned int>(dof_indices.size());
256 348 : const unsigned int n_qp = qrule->n_points();
257 :
258 54342 : for (unsigned int qp = 0; qp < n_qp; qp++)
259 : {
260 58548 : std::vector<Real> psi = legepoly(dim, element_order, q_point[qp], matsize);
261 8372 : const unsigned int psi_size = cast_int<unsigned int>(psi.size());
262 :
263 345492 : for (unsigned int i = 0; i < matsize; i++)
264 2065960 : for (unsigned int j = 0; j < matsize; j++)
265 2361316 : Kp(i, j) += JxW[qp] * psi[i] * psi[j];
266 :
267 4186 : Number u_h = libMesh::zero;
268 466452 : for (unsigned int i = 0; i < n_dofs; i++)
269 520412 : u_h += (*phi)[i][qp] * system.current_solution(dof_indices[i]);
270 :
271 345492 : for (unsigned int i = 0; i < psi_size; i++)
272 369212 : F(i) += JxW[qp] * u_h * psi[i];
273 : }
274 :
275 4166 : Kp.lu_solve(F, Pu_h);
276 :
277 : // Generate index to total degree map. Total_degree_per_index[i] gives the degree ith basis function
278 696 : std::vector<unsigned int> total_degree_per_index;
279 :
280 19936 : for (unsigned int poly_deg = 0; poly_deg <= static_cast<unsigned int>(element_order); ++poly_deg)
281 : {
282 15770 : switch (dim)
283 : {
284 0 : case 3:
285 0 : for (int i = poly_deg; i >= 0; --i)
286 0 : for (int j = poly_deg - i; j >= 0; --j)
287 : {
288 0 : int k = poly_deg - i - j;
289 0 : total_degree_per_index.push_back(i + j + k);
290 0 : }
291 0 : break;
292 :
293 7560 : case 2:
294 22680 : for (int i = poly_deg; i >= 0; --i)
295 : {
296 2520 : int j = poly_deg - i;
297 15120 : total_degree_per_index.push_back(i + j);
298 1260 : }
299 630 : break;
300 :
301 8210 : case 1:
302 8210 : total_degree_per_index.push_back(poly_deg);
303 688 : break;
304 :
305 0 : default:
306 0 : libmesh_error_msg("Invalid dimension dim " << dim);
307 : }
308 : }
309 :
310 : // Group coefficients |c_k| by degree k
311 696 : std::map<unsigned int, std::vector<Number>> coeff_by_degree;
312 :
313 27496 : for (unsigned int i = 0; i < Pu_h.size(); ++i)
314 : {
315 23330 : unsigned int degree = total_degree_per_index[i];
316 : // Excluding the constant mode i.e, zeroth order coefficient as we use ln|order| later
317 : // Constant order term often dominates the scale and doesn't inform about regularity.
318 23330 : if (degree > 0)
319 19164 : coeff_by_degree[degree].push_back(Pu_h(i));
320 : }
321 :
322 : // Compute L2 norm of each group |c_k|
323 696 : std::vector<unsigned int> degrees;
324 696 : std::vector<double> norms;
325 :
326 15770 : for (const auto & pair : coeff_by_degree)
327 : {
328 11604 : unsigned int deg = pair.first;
329 970 : const std::vector<Number> & coeffs = pair.second;
330 :
331 11604 : double norm = 0.;
332 30768 : for (const auto & c : coeffs)
333 20764 : norm += std::norm(c);
334 :
335 11604 : norm = std::sqrt(norm);
336 :
337 11604 : degrees.push_back(deg);
338 11604 : norms.push_back(norm);
339 : }
340 :
341 :
342 : // Collect log |c_k| and log |k|
343 696 : std::vector<double> log_norms, log_degrees;
344 :
345 16118 : for (unsigned int i = 0; i < degrees.size(); ++i)
346 : {
347 : // filter tiny terms
348 12574 : if (norms[i] > 1e-12)
349 : {
350 11604 : log_degrees.push_back(std::log(static_cast<double>(degrees[i])));
351 12574 : log_norms.push_back(std::log(norms[i]));
352 : }
353 : }
354 :
355 :
356 : // Fit log-log line - we use least-squares fit
357 348 : double Sx = 0., Sy = 0., Sxx = 0., Sxy = 0.;
358 696 : const size_t N = log_degrees.size();
359 :
360 15770 : for (size_t i = 0; i < N; ++i)
361 : {
362 11604 : Sx += log_degrees[i];
363 11604 : Sy += log_norms[i];
364 11604 : Sxx += log_degrees[i] * log_degrees[i];
365 11604 : Sxy += log_degrees[i] * log_norms[i];
366 : }
367 :
368 4166 : const double regularity = -compute_slope(N, Sx, Sy, Sxx, Sxy);
369 : // const double intercept = (Sy - slope * Sx) / N;
370 :
371 348 : Threads::spin_mutex::scoped_lock acquire(Threads::spin_mtx);
372 4514 : smoothness_per_cell[e_id] = regularity;
373 6940 : } // end variables loop
374 :
375 : } // end element loop
376 :
377 2414 : } // End () operator definition
378 :
379 : } // namespace libMesh
|
