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