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 : // 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 748 : 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 780 : for (auto & child : elem->child_ref_range())
57 552 : this->add_projection(system, &child, var);
58 228 : return;
59 : }
60 :
61 : // The DofMap for this system
62 132 : const DofMap & dof_map = system.get_dof_map();
63 :
64 520 : const FEContinuity cont = fe->get_continuity();
65 :
66 520 : fe->reinit(elem);
67 :
68 520 : 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 520 : FEMap::inverse_map (system.get_mesh().mesh_dimension(), coarse,
74 520 : *xyz_values, coarse_qpoints);
75 :
76 520 : fe_coarse->reinit(coarse, &coarse_qpoints);
77 :
78 : const unsigned int n_coarse_dofs =
79 520 : cast_int<unsigned int>(phi_coarse->size());
80 :
81 520 : if (Uc.size() == 0)
82 : {
83 146 : Ke.resize(n_coarse_dofs, n_coarse_dofs);
84 50 : Ke.zero();
85 146 : Fe.resize(n_coarse_dofs);
86 50 : Fe.zero();
87 146 : 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 5602 : 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 43930 : for (unsigned int i=0; i != n_dofs; i++)
101 : {
102 38848 : dof_id_type dof_num = dof_indices[i];
103 48594 : val += (*phi)[i][qp] *
104 38848 : system.current_solution(dof_num);
105 38848 : if (cont == C_ZERO || cont == C_ONE)
106 38848 : grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
107 38848 : if (cont == C_ONE)
108 10744 : hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
109 : }
110 :
111 : // The projection matrix and vector
112 44170 : for (auto i : index_range(Fe))
113 : {
114 39088 : Fe(i) += (*JxW)[qp] *
115 58700 : (*phi_coarse)[i][qp]*val;
116 39088 : if (cont == C_ZERO || cont == C_ONE)
117 39272 : Fe(i) += (*JxW)[qp] *
118 48894 : (grad*(*dphi_coarse)[i][qp]);
119 39088 : if (cont == C_ONE)
120 11144 : Fe(i) += (*JxW)[qp] *
121 13758 : hess.contract((*d2phi_coarse)[i][qp]);
122 :
123 359160 : for (auto j : index_range(Fe))
124 : {
125 400196 : Ke(i,j) += (*JxW)[qp] *
126 480320 : (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
127 320072 : if (cont == C_ZERO || cont == C_ONE)
128 240780 : Ke(i,j) += (*JxW)[qp] *
129 560444 : (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
130 320072 : if (cont == C_ONE)
131 70632 : Ke(i,j) += (*JxW)[qp] *
132 121780 : ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
133 : }
134 : }
135 : }
136 : }
137 :
138 119 : 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 119 : const unsigned int dim = mesh.mesh_dimension();
147 :
148 : // The number of variables in the system
149 119 : 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 119 : 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 119 : 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 153 : SmoothnessEstimator estimate_smoothness;
177 119 : 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 153 : std::vector<ErrorVectorReal> h_error_per_cell(mesh.max_elem_id(), 0.);
182 119 : std::vector<ErrorVectorReal> p_error_per_cell(mesh.max_elem_id(), 0.);
183 :
184 : // Loop over all the variables in the system
185 238 : for (unsigned int var=0; var<n_vars; var++)
186 : {
187 : // Possibly skip this variable
188 119 : if (!component_scale.empty())
189 153 : 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 170 : fe = FEBase::build (dim, fe_type);
197 170 : fe_coarse = FEBase::build (dim, fe_type);
198 :
199 : // Any cached coarse element results have expired
200 119 : coarse = nullptr;
201 34 : unsigned int cached_coarse_p_level = 0;
202 :
203 119 : 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 204 : 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 153 : 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 119 : JxW = &(fe->get_JxW());
217 119 : xyz_values = &(fe->get_xyz());
218 :
219 : // The shape functions
220 119 : phi = &(fe->get_phi());
221 119 : phi_coarse = &(fe_coarse->get_phi());
222 :
223 : // The shape function derivatives
224 119 : if (cont == C_ZERO || cont == C_ONE)
225 : {
226 119 : dphi = &(fe->get_dphi());
227 119 : dphi_coarse = &(fe_coarse->get_dphi());
228 : }
229 :
230 : // The shape function second derivatives
231 119 : if (cont == C_ONE)
232 : {
233 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
234 49 : d2phi = &(fe->get_d2phi());
235 49 : 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 1248 : 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 696 : if (elem->refinement_flag() != Elem::REFINE)
248 291 : 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 457 : if (elem->parent() &&
255 288 : (coarse != elem->parent() ||
256 44 : cached_coarse_p_level != elem->p_level()))
257 : {
258 146 : Uc.resize(0);
259 :
260 196 : coarse = elem->parent();
261 196 : 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 196 : this->add_projection(system, coarse, var);
267 :
268 196 : coarse->hack_p_level(old_parent_level);
269 :
270 : // Solve the h-coarsening projection problem
271 196 : Ke.cholesky_solve(Fe, Uc);
272 : }
273 :
274 308 : fe->reinit(elem);
275 :
276 : // Get the DOF indices for the fine element
277 308 : 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 308 : 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 385 : if (elem->p_level() == 0)
298 : {
299 43 : unsigned int n_vertices = 0;
300 1048 : for (unsigned int n = 0; n != n_nodes; ++n)
301 876 : if (elem->is_vertex(n))
302 : {
303 464 : n_vertices++;
304 464 : const Node & node = elem->node_ref(n);
305 348 : average_val += system.current_solution
306 464 : (node.dof_number(sys_num,var,0));
307 : }
308 172 : 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 170 : fe_coarse->reinit(elem, &(qrule->get_points()));
316 :
317 : const unsigned int n_coarse_dofs =
318 136 : cast_int<unsigned int>(phi_coarse->size());
319 :
320 136 : elem->hack_p_level(old_elem_level);
321 :
322 136 : Ke.resize(n_coarse_dofs, n_coarse_dofs);
323 34 : Ke.zero();
324 102 : Fe.resize(n_coarse_dofs);
325 34 : Fe.zero();
326 :
327 : // Loop over the quadrature points
328 1126 : 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 7446 : for (unsigned int i=0; i != n_dofs; i++)
336 : {
337 6456 : dof_id_type dof_num = dof_indices[i];
338 8064 : val += (*phi)[i][qp] *
339 6456 : system.current_solution(dof_num);
340 6456 : if (cont == C_ZERO || cont == C_ONE)
341 6456 : grad.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
342 6456 : if (cont == C_ONE)
343 6456 : hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
344 : }
345 :
346 : // The projection matrix and vector
347 6456 : for (auto i : index_range(Fe))
348 : {
349 5466 : Fe(i) += (*JxW)[qp] *
350 8188 : (*phi_coarse)[i][qp]*val;
351 5466 : if (cont == C_ZERO || cont == C_ONE)
352 5466 : Fe(i) += (*JxW)[qp] *
353 8188 : grad * (*dphi_coarse)[i][qp];
354 5466 : if (cont == C_ONE)
355 5466 : Fe(i) += (*JxW)[qp] *
356 6827 : hess.contract((*d2phi_coarse)[i][qp]);
357 :
358 37548 : for (auto j : index_range(Fe))
359 : {
360 40063 : Ke(i,j) += (*JxW)[qp] *
361 48044 : (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
362 32082 : if (cont == C_ZERO || cont == C_ONE)
363 24101 : Ke(i,j) += (*JxW)[qp] *
364 56025 : (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
365 32082 : if (cont == C_ONE)
366 32082 : Ke(i,j) += (*JxW)[qp] *
367 56025 : ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
368 : }
369 : }
370 : }
371 :
372 : // Solve the p-coarsening projection problem
373 136 : Ke.cholesky_solve(Fe, Up);
374 : }
375 :
376 : // loop over the integration points on the fine element
377 2782 : 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 19522 : for (unsigned int i=0; i<n_dofs; i++)
383 : {
384 17048 : const dof_id_type dof_num = dof_indices[i];
385 21304 : value_error += (*phi)[i][qp] *
386 17048 : system.current_solution(dof_num);
387 17048 : if (cont == C_ZERO || cont == C_ONE)
388 17048 : grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
389 17048 : if (cont == C_ONE)
390 7976 : hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
391 : }
392 2474 : if (elem->p_level() == 0)
393 : {
394 1113 : value_error -= average_val;
395 : }
396 : else
397 : {
398 6456 : for (auto i : index_range(Up))
399 : {
400 8188 : value_error -= (*phi_coarse)[i][qp] * Up(i);
401 5466 : if (cont == C_ZERO || cont == C_ONE)
402 6827 : grad_error.subtract_scaled((*dphi_coarse)[i][qp], Up(i));
403 5466 : if (cont == C_ONE)
404 6827 : hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Up(i));
405 : }
406 : }
407 :
408 3710 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
409 3092 : (component_scale[var] *
410 3092 : (*JxW)[qp] * TensorTools::norm_sq(value_error));
411 2474 : if (cont == C_ZERO || cont == C_ONE)
412 2474 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
413 2474 : (component_scale[var] *
414 2474 : (*JxW)[qp] * grad_error.norm_sq());
415 2474 : if (cont == C_ONE)
416 1370 : p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
417 1370 : (component_scale[var] *
418 1370 : (*JxW)[qp] * hessian_error.norm_sq());
419 : }
420 :
421 : // Calculate this variable's contribution to the h
422 : // refinement error
423 :
424 385 : if (!elem->parent())
425 : {
426 : // For now, we'll always start with an h refinement
427 25 : h_error_per_cell[e_id] =
428 5 : std::numeric_limits<ErrorVectorReal>::max() / 2;
429 : }
430 : else
431 : {
432 288 : FEMap::inverse_map (dim, coarse, *xyz_values,
433 288 : coarse_qpoints);
434 :
435 288 : unsigned int old_parent_level = coarse->p_level();
436 288 : coarse->hack_p_level(elem->p_level());
437 :
438 288 : fe_coarse->reinit(coarse, &coarse_qpoints);
439 :
440 288 : coarse->hack_p_level(old_parent_level);
441 :
442 : // The number of DOFS on the coarse element
443 : unsigned int n_coarse_dofs =
444 288 : cast_int<unsigned int>(phi_coarse->size());
445 :
446 : // Loop over the quadrature points
447 2534 : 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 17438 : for (unsigned int i=0; i != n_dofs; ++i)
455 : {
456 15192 : const dof_id_type dof_num = dof_indices[i];
457 18984 : value_error += (*phi)[i][qp] *
458 15192 : system.current_solution(dof_num);
459 15192 : if (cont == C_ZERO || cont == C_ONE)
460 15192 : grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
461 15192 : if (cont == C_ONE)
462 7896 : hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
463 : }
464 :
465 17438 : for (unsigned int i=0; i != n_coarse_dofs; ++i)
466 : {
467 22776 : value_error -= (*phi_coarse)[i][qp] * Uc(i);
468 15192 : if (cont == C_ZERO || cont == C_ONE)
469 18984 : grad_error.subtract_scaled((*dphi_coarse)[i][qp], Uc(i));
470 15192 : if (cont == C_ONE)
471 9864 : hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Uc(i));
472 : }
473 :
474 3368 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
475 2807 : (component_scale[var] *
476 2807 : (*JxW)[qp] * TensorTools::norm_sq(value_error));
477 2246 : if (cont == C_ZERO || cont == C_ONE)
478 2246 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
479 2246 : (component_scale[var] *
480 2246 : (*JxW)[qp] * grad_error.norm_sq());
481 2246 : if (cont == C_ONE)
482 1350 : h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
483 1350 : (component_scale[var] *
484 1350 : (*JxW)[qp] * hessian_error.norm_sq());
485 : }
486 :
487 : }
488 51 : }
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 1248 : 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 696 : if (elem->refinement_flag() != Elem::REFINE)
501 291 : 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 616 : 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 308 : dofs_per_elem +=
518 308 : FEInterface::n_dofs(fe_type, elem);
519 :
520 : // Compute number of DOFs for elem at current p_level() + 1
521 308 : dofs_per_p_elem +=
522 385 : FEInterface::n_dofs(fe_type, elem->p_level() + 1, elem);
523 : }
524 :
525 : const unsigned int new_h_dofs = dofs_per_elem *
526 308 : (elem->n_children() - 1);
527 :
528 308 : 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 385 : std::sqrt(p_error_per_cell[e_id]) * p_weight / new_p_dofs;
541 : const Real h_value =
542 385 : std::sqrt(h_error_per_cell[e_id]) /
543 308 : static_cast<Real>(new_h_dofs);
544 385 : if (smoothness[elem->id()] > smoothness.mean() && p_value > h_value)
545 : {
546 111 : elem->set_p_refinement_flag(Elem::REFINE);
547 54 : elem->set_refinement_flag(Elem::DO_NOTHING);
548 : }
549 51 : }
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 119 : MeshRefinement(mesh).make_flags_parallel_consistent();
557 119 : }
558 :
559 : } // namespace libMesh
560 :
561 : #endif // #ifdef LIBMESH_ENABLE_AMR
|
