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 : // C++ includes
19 : #include <limits> // for std::numeric_limits::max
20 : #include <math.h> // for sqrt
21 :
22 :
23 : // Local Includes
24 : #include "libmesh/hp_coarsentest.h"
25 : #include "libmesh/dense_matrix.h"
26 : #include "libmesh/dense_vector.h"
27 : #include "libmesh/dof_map.h"
28 : #include "libmesh/fe_base.h"
29 : #include "libmesh/fe_interface.h"
30 : #include "libmesh/libmesh_logging.h"
31 : #include "libmesh/elem.h"
32 : #include "libmesh/error_vector.h"
33 : #include "libmesh/mesh_base.h"
34 : #include "libmesh/mesh_refinement.h"
35 : #include "libmesh/quadrature.h"
36 : #include "libmesh/system.h"
37 : #include "libmesh/tensor_value.h"
38 : #include "libmesh/smoothness_estimator.h"
39 :
40 : #ifdef LIBMESH_ENABLE_AMR
41 :
42 : namespace libMesh
43 : {
44 :
45 : //-----------------------------------------------------------------
46 : // HPCoarsenTest implementations
47 :
48 2509 : void HPCoarsenTest::add_projection(const System & system,
49 : const Elem * elem,
50 : unsigned int var)
51 : {
52 : // If we have children, we need to add their projections instead
53 190 : if (!elem->active())
54 : {
55 58 : libmesh_assert(!elem->subactive());
56 2604 : for (auto & child : elem->child_ref_range())
57 1832 : this->add_projection(system, &child, var);
58 772 : return;
59 : }
60 :
61 : // The DofMap for this system
62 132 : const DofMap & dof_map = system.get_dof_map();
63 :
64 1737 : const FEContinuity cont = fe->get_continuity();
65 :
66 1737 : fe->reinit(elem);
67 :
68 1737 : dof_map.dof_indices(elem, dof_indices, var);
69 :
70 : const unsigned int n_dofs =
71 264 : cast_int<unsigned int>(dof_indices.size());
72 :
73 1737 : FEMap::inverse_map (system.get_mesh().mesh_dimension(), coarse,
74 1737 : *xyz_values, coarse_qpoints);
75 :
76 1737 : fe_coarse->reinit(coarse, &coarse_qpoints);
77 :
78 : const unsigned int n_coarse_dofs =
79 1737 : cast_int<unsigned int>(phi_coarse->size());
80 :
81 1737 : if (Uc.size() == 0)
82 : {
83 627 : Ke.resize(n_coarse_dofs, n_coarse_dofs);
84 50 : Ke.zero();
85 627 : Fe.resize(n_coarse_dofs);
86 50 : Fe.zero();
87 627 : Uc.resize(n_coarse_dofs);
88 50 : Uc.zero();
89 : }
90 132 : libmesh_assert_equal_to (Uc.size(), phi_coarse->size());
91 :
92 : // Loop over the quadrature points
93 17954 : for (auto qp : make_range(qrule->n_points()))
94 : {
95 : // The solution value at the quadrature point
96 1280 : Number val = libMesh::zero;
97 1280 : Gradient grad;
98 2560 : Tensor hess;
99 :
100 137477 : for (unsigned int i=0; i != n_dofs; i++)
101 : {
102 121260 : dof_id_type dof_num = dof_indices[i];
103 131006 : val += (*phi)[i][qp] *
104 121260 : system.current_solution(dof_num);
105 121260 : if (cont == C_ZERO || cont == C_ONE)
106 121260 : grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
107 121260 : if (cont == C_ONE)
108 36444 : hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
109 : }
110 :
111 : // The projection matrix and vector
112 138220 : for (auto i : index_range(Fe))
113 : {
114 122003 : Fe(i) += (*JxW)[qp] *
115 141615 : (*phi_coarse)[i][qp]*val;
116 122003 : if (cont == C_ZERO || cont == C_ONE)
117 122187 : Fe(i) += (*JxW)[qp] *
118 131809 : (grad*(*dphi_coarse)[i][qp]);
119 122003 : if (cont == C_ONE)
120 37347 : Fe(i) += (*JxW)[qp] *
121 39961 : hess.contract((*d2phi_coarse)[i][qp]);
122 :
123 1107950 : for (auto j : index_range(Fe))
124 : {
125 1066071 : Ke(i,j) += (*JxW)[qp] *
126 1146195 : (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
127 985947 : if (cont == C_ZERO || cont == C_ONE)
128 906655 : Ke(i,j) += (*JxW)[qp] *
129 1226319 : (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
130 985947 : if (cont == C_ONE)
131 234595 : Ke(i,j) += (*JxW)[qp] *
132 285743 : ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
133 : }
134 : }
135 : }
136 : }
137 :
138 1207 : void HPCoarsenTest::select_refinement (System & system)
139 : {
140 68 : LOG_SCOPE("select_refinement()", "HPCoarsenTest");
141 :
142 : // The current mesh
143 68 : MeshBase & mesh = system.get_mesh();
144 :
145 : // The dimensionality of the mesh
146 1207 : const unsigned int dim = mesh.mesh_dimension();
147 :
148 : // The number of variables in the system
149 34 : const unsigned int n_vars = system.n_vars();
150 :
151 : // The DofMap for this system
152 34 : const DofMap & dof_map = system.get_dof_map();
153 :
154 : // The system number (for doing bad hackery)
155 68 : const unsigned int sys_num = system.number();
156 :
157 : // Check for a valid component_scale
158 1207 : if (!component_scale.empty())
159 : {
160 0 : libmesh_error_msg_if(component_scale.size() != n_vars,
161 : "ERROR: component_scale is the wrong size:\n"
162 : << " component_scale.size()="
163 : << component_scale.size()
164 : << "\n n_vars="
165 : << n_vars);
166 : }
167 : else
168 : {
169 : // No specified scaling. Scale all variables by one.
170 1207 : component_scale.resize (n_vars, 1.0);
171 : }
172 :
173 : // Estimates smoothness of solution on each cell
174 : // via Legendre coefficient decay rate.
175 68 : ErrorVector smoothness;
176 1241 : SmoothnessEstimator estimate_smoothness;
177 1207 : estimate_smoothness.estimate_smoothness(system, smoothness);
178 :
179 : // Resize the error_per_cell vectors to handle
180 : // the number of elements, initialize them to 0.
181 1241 : std::vector<ErrorVectorReal> h_error_per_cell(mesh.max_elem_id(), 0.);
182 1207 : std::vector<ErrorVectorReal> p_error_per_cell(mesh.max_elem_id(), 0.);
183 :
184 : // Loop over all the variables in the system
185 2414 : for (unsigned int var=0; var<n_vars; var++)
186 : {
187 : // Possibly skip this variable
188 1207 : if (!component_scale.empty())
189 1241 : if (component_scale[var] == 0.0) continue;
190 :
191 : // The type of finite element to use for this variable
192 34 : const FEType & fe_type = dof_map.variable_type (var);
193 :
194 : // Finite element objects for a fine (and probably a coarse)
195 : // element will be needed
196 2346 : fe = FEBase::build (dim, fe_type);
197 2346 : fe_coarse = FEBase::build (dim, fe_type);
198 :
199 : // Any cached coarse element results have expired
200 1207 : coarse = nullptr;
201 34 : unsigned int cached_coarse_p_level = 0;
202 :
203 1207 : const FEContinuity cont = fe->get_continuity();
204 34 : libmesh_assert (cont == DISCONTINUOUS || cont == C_ZERO ||
205 : cont == C_ONE);
206 :
207 : // Build an appropriate quadrature rule
208 2380 : qrule = fe_type.default_quadrature_rule(dim, _extra_order);
209 :
210 : // Tell the refined finite element about the quadrature
211 : // rule. The coarse finite element need not know about it
212 1241 : fe->attach_quadrature_rule (qrule.get());
213 :
214 : // We will always do the integration
215 : // on the fine elements. Get their Jacobian values, etc..
216 1207 : JxW = &(fe->get_JxW());
217 1207 : xyz_values = &(fe->get_xyz());
218 :
219 : // The shape functions
220 1207 : phi = &(fe->get_phi());
221 1207 : phi_coarse = &(fe_coarse->get_phi());
222 :
223 : // The shape function derivatives
224 1207 : if (cont == C_ZERO || cont == C_ONE)
225 : {
226 1207 : dphi = &(fe->get_dphi());
227 1207 : dphi_coarse = &(fe_coarse->get_dphi());
228 : }
229 :
230 : // The shape function second derivatives
231 1207 : if (cont == C_ONE)
232 : {
233 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
234 497 : d2phi = &(fe->get_d2phi());
235 497 : d2phi_coarse = &(fe_coarse->get_d2phi());
236 : #else
237 : libmesh_error_msg("Minimization of H2 error without second derivatives is not possible.");
238 : #endif
239 : }
240 :
241 : // Iterate over all the active elements in the mesh
242 : // that live on this processor.
243 6198 : for (auto & elem : mesh.active_local_element_ptr_range())
244 : {
245 : // We're only checking elements that are already flagged for h
246 : // refinement
247 2083 : if (elem->refinement_flag() != Elem::REFINE)
248 1067 : continue;
249 :
250 154 : const dof_id_type e_id = elem->id();
251 :
252 : // Find the projection onto the parent element,
253 : // if necessary
254 1068 : if (elem->parent() &&
255 859 : (coarse != elem->parent() ||
256 44 : cached_coarse_p_level != elem->p_level()))
257 : {
258 627 : Uc.resize(0);
259 :
260 677 : coarse = elem->parent();
261 677 : cached_coarse_p_level = elem->p_level();
262 :
263 100 : unsigned int old_parent_level = coarse->p_level();
264 50 : coarse->hack_p_level(elem->p_level());
265 :
266 677 : this->add_projection(system, coarse, var);
267 :
268 677 : coarse->hack_p_level(old_parent_level);
269 :
270 : // Solve the h-coarsening projection problem
271 677 : Ke.cholesky_solve(Fe, Uc);
272 : }
273 :
274 919 : fe->reinit(elem);
275 :
276 : // Get the DOF indices for the fine element
277 919 : dof_map.dof_indices (elem, dof_indices, var);
278 :
279 : // The number of quadrature points
280 77 : const unsigned int n_qp = qrule->n_points();
281 :
282 : // The number of DOFS on the fine element
283 : const unsigned int n_dofs =
284 154 : cast_int<unsigned int>(dof_indices.size());
285 :
286 : // The number of nodes on the fine element
287 919 : const unsigned int n_nodes = elem->n_nodes();
288 :
289 : // The average element value (used as an ugly hack
290 : // when we have nothing p-coarsened to compare to)
291 : // Real average_val = 0.;
292 154 : Number average_val = 0.;
293 :
294 : // Calculate this variable's contribution to the p
295 : // refinement error
296 :
297 996 : if (elem->p_level() == 0)
298 : {
299 43 : unsigned int n_vertices = 0;
300 3136 : for (unsigned int n = 0; n != n_nodes; ++n)
301 2622 : if (elem->is_vertex(n))
302 : {
303 1388 : n_vertices++;
304 1388 : const Node & node = elem->node_ref(n);
305 1272 : average_val += system.current_solution
306 1388 : (node.dof_number(sys_num,var,0));
307 : }
308 514 : average_val /= n_vertices;
309 : }
310 : else
311 : {
312 34 : unsigned int old_elem_level = elem->p_level();
313 34 : elem->hack_p_level(old_elem_level - 1);
314 :
315 439 : fe_coarse->reinit(elem, &(qrule->get_points()));
316 :
317 : const unsigned int n_coarse_dofs =
318 405 : cast_int<unsigned int>(phi_coarse->size());
319 :
320 405 : elem->hack_p_level(old_elem_level);
321 :
322 405 : Ke.resize(n_coarse_dofs, n_coarse_dofs);
323 34 : Ke.zero();
324 371 : Fe.resize(n_coarse_dofs);
325 34 : Fe.zero();
326 :
327 : // Loop over the quadrature points
328 3349 : for (auto qp : make_range(qrule->n_points()))
329 : {
330 : // The solution value at the quadrature point
331 247 : Number val = libMesh::zero;
332 247 : Gradient grad;
333 494 : Tensor hess;
334 :
335 22122 : for (unsigned int i=0; i != n_dofs; i++)
336 : {
337 19178 : dof_id_type dof_num = dof_indices[i];
338 20786 : val += (*phi)[i][qp] *
339 19178 : system.current_solution(dof_num);
340 19178 : if (cont == C_ZERO || cont == C_ONE)
341 19178 : grad.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
342 19178 : if (cont == C_ONE)
343 19178 : hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
344 : }
345 :
346 : // The projection matrix and vector
347 19178 : for (auto i : index_range(Fe))
348 : {
349 16234 : Fe(i) += (*JxW)[qp] *
350 18956 : (*phi_coarse)[i][qp]*val;
351 16234 : if (cont == C_ZERO || cont == C_ONE)
352 16234 : Fe(i) += (*JxW)[qp] *
353 18956 : grad * (*dphi_coarse)[i][qp];
354 16234 : if (cont == C_ONE)
355 16234 : Fe(i) += (*JxW)[qp] *
356 17595 : hess.contract((*d2phi_coarse)[i][qp]);
357 :
358 111492 : for (auto j : index_range(Fe))
359 : {
360 103239 : Ke(i,j) += (*JxW)[qp] *
361 111220 : (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
362 95258 : if (cont == C_ZERO || cont == C_ONE)
363 87277 : Ke(i,j) += (*JxW)[qp] *
364 119201 : (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
365 95258 : if (cont == C_ONE)
366 95258 : Ke(i,j) += (*JxW)[qp] *
367 119201 : ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
368 : }
369 : }
370 : }
371 :
372 : // Solve the p-coarsening projection problem
373 405 : Ke.cholesky_solve(Fe, Up);
374 : }
375 :
376 : // loop over the integration points on the fine element
377 8305 : for (unsigned int qp=0; qp<n_qp; qp++)
378 : {
379 618 : Number value_error = 0.;
380 618 : Gradient grad_error;
381 1236 : Tensor hessian_error;
382 58300 : for (unsigned int i=0; i<n_dofs; i++)
383 : {
384 50914 : const dof_id_type dof_num = dof_indices[i];
385 55170 : value_error += (*phi)[i][qp] *
386 50914 : system.current_solution(dof_num);
387 50914 : if (cont == C_ZERO || cont == C_ONE)
388 50914 : grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
389 50914 : if (cont == C_ONE)
390 23698 : hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
391 : }
392 7386 : if (elem->p_level() == 0)
393 : {
394 4071 : value_error -= average_val;
395 : }
396 : else
397 : {
398 19178 : for (auto i : index_range(Up))
399 : {
400 18956 : value_error -= (*phi_coarse)[i][qp] * Up(i);
401 16234 : if (cont == C_ZERO || cont == C_ONE)
402 17595 : grad_error.subtract_scaled((*dphi_coarse)[i][qp], Up(i));
403 16234 : if (cont == C_ONE)
404 17595 : hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Up(i));
405 : }
406 : }
407 :
408 8622 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
409 8004 : (component_scale[var] *
410 8004 : (*JxW)[qp] * TensorTools::norm_sq(value_error));
411 7386 : if (cont == C_ZERO || cont == C_ONE)
412 7386 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
413 7386 : (component_scale[var] *
414 7386 : (*JxW)[qp] * grad_error.norm_sq());
415 7386 : if (cont == C_ONE)
416 4074 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
417 4074 : (component_scale[var] *
418 4074 : (*JxW)[qp] * hessian_error.norm_sq());
419 : }
420 :
421 : // Calculate this variable's contribution to the h
422 : // refinement error
423 :
424 996 : if (!elem->parent())
425 : {
426 : // For now, we'll always start with an h refinement
427 65 : h_error_per_cell[e_id] =
428 5 : std::numeric_limits<ErrorVectorReal>::max() / 2;
429 : }
430 : else
431 : {
432 859 : FEMap::inverse_map (dim, coarse, *xyz_values,
433 859 : coarse_qpoints);
434 :
435 859 : unsigned int old_parent_level = coarse->p_level();
436 859 : coarse->hack_p_level(elem->p_level());
437 :
438 859 : fe_coarse->reinit(coarse, &coarse_qpoints);
439 :
440 859 : coarse->hack_p_level(old_parent_level);
441 :
442 : // The number of DOFS on the coarse element
443 : unsigned int n_coarse_dofs =
444 859 : cast_int<unsigned int>(phi_coarse->size());
445 :
446 : // Loop over the quadrature points
447 7561 : for (unsigned int qp=0; qp<n_qp; qp++)
448 : {
449 : // The solution difference at the quadrature point
450 561 : Number value_error = libMesh::zero;
451 561 : Gradient grad_error;
452 1122 : Tensor hessian_error;
453 :
454 52048 : for (unsigned int i=0; i != n_dofs; ++i)
455 : {
456 45346 : const dof_id_type dof_num = dof_indices[i];
457 49138 : value_error += (*phi)[i][qp] *
458 45346 : system.current_solution(dof_num);
459 45346 : if (cont == C_ZERO || cont == C_ONE)
460 45346 : grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
461 45346 : if (cont == C_ONE)
462 23458 : hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
463 : }
464 :
465 52048 : for (unsigned int i=0; i != n_coarse_dofs; ++i)
466 : {
467 52930 : value_error -= (*phi_coarse)[i][qp] * Uc(i);
468 45346 : if (cont == C_ZERO || cont == C_ONE)
469 49138 : grad_error.subtract_scaled((*dphi_coarse)[i][qp], Uc(i));
470 45346 : if (cont == C_ONE)
471 25426 : hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Uc(i));
472 : }
473 :
474 7824 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
475 7263 : (component_scale[var] *
476 7263 : (*JxW)[qp] * TensorTools::norm_sq(value_error));
477 6702 : if (cont == C_ZERO || cont == C_ONE)
478 6702 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
479 6702 : (component_scale[var] *
480 6702 : (*JxW)[qp] * grad_error.norm_sq());
481 6702 : if (cont == C_ONE)
482 4014 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
483 4014 : (component_scale[var] *
484 4014 : (*JxW)[qp] * hessian_error.norm_sq());
485 : }
486 :
487 : }
488 1139 : }
489 : }
490 :
491 : // Now that we've got our approximations for p_error and h_error, let's see
492 : // if we want to switch any h refinement flags to p refinement
493 :
494 : // Iterate over all the active elements in the mesh
495 : // that live on this processor.
496 6198 : for (auto & elem : mesh.active_local_element_ptr_range())
497 : {
498 : // We're only checking elements that are already flagged for h
499 : // refinement
500 2083 : if (elem->refinement_flag() != Elem::REFINE)
501 1067 : continue;
502 :
503 154 : const dof_id_type e_id = elem->id();
504 :
505 77 : unsigned int dofs_per_elem = 0, dofs_per_p_elem = 0;
506 :
507 : // Loop over all the variables in the system
508 1838 : for (unsigned int var=0; var<n_vars; var++)
509 : {
510 : // The type of finite element to use for this variable
511 77 : const FEType & fe_type = dof_map.variable_type (var);
512 :
513 : // FIXME: we're overestimating the number of DOFs added by h
514 : // refinement
515 :
516 : // Compute number of DOFs for elem at current p_level()
517 919 : dofs_per_elem +=
518 919 : FEInterface::n_dofs(fe_type, elem);
519 :
520 : // Compute number of DOFs for elem at current p_level() + 1
521 919 : dofs_per_p_elem +=
522 996 : FEInterface::n_dofs(fe_type, elem->p_level() + 1, elem);
523 : }
524 :
525 : const unsigned int new_h_dofs = dofs_per_elem *
526 919 : (elem->n_children() - 1);
527 :
528 919 : const unsigned int new_p_dofs = dofs_per_p_elem -
529 : dofs_per_elem;
530 :
531 : /*
532 : libMesh::err << "Cell " << e_id << ": h = " << elem->hmax()
533 : << ", p = " << elem->p_level() + 1 << "," << std::endl
534 : << " h_error = " << h_error_per_cell[e_id]
535 : << ", p_error = " << p_error_per_cell[e_id] << std::endl
536 : << " new_h_dofs = " << new_h_dofs
537 : << ", new_p_dofs = " << new_p_dofs << std::endl;
538 : */
539 : const Real p_value =
540 996 : std::sqrt(p_error_per_cell[e_id]) * p_weight / new_p_dofs;
541 : const Real h_value =
542 996 : std::sqrt(h_error_per_cell[e_id]) /
543 919 : static_cast<Real>(new_h_dofs);
544 996 : if (smoothness[elem->id()] > smoothness.mean() && p_value > h_value)
545 : {
546 335 : elem->set_p_refinement_flag(Elem::REFINE);
547 56 : elem->set_refinement_flag(Elem::DO_NOTHING);
548 : }
549 1139 : }
550 :
551 : // libMesh::MeshRefinement will now assume that users have set
552 : // refinement flags consistently on all processors, but all we've
553 : // done so far is set refinement flags on local elements. Let's
554 : // make sure that flags on geometrically ghosted elements are all
555 : // set to whatever their owners decided.
556 1207 : MeshRefinement(mesh).make_flags_parallel_consistent();
557 1207 : }
558 :
559 : } // namespace libMesh
560 :
561 : #endif // #ifdef LIBMESH_ENABLE_AMR
|
