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 : // Local Includes
20 : #include "libmesh/libmesh_common.h"
21 : #include "libmesh/jump_error_estimator.h"
22 : #include "libmesh/dof_map.h"
23 : #include "libmesh/elem.h"
24 : #include "libmesh/error_vector.h"
25 : #include "libmesh/fe_base.h"
26 : #include "libmesh/fe_interface.h"
27 : #include "libmesh/fem_context.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/quadrature_gauss.h"
31 : #include "libmesh/system.h"
32 : #include "libmesh/dense_vector.h"
33 : #include "libmesh/numeric_vector.h"
34 : #include "libmesh/int_range.h"
35 :
36 : // C++ Includes
37 : #include <algorithm> // for std::fill
38 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
39 : #include <cmath> // for sqrt
40 : #include <memory>
41 :
42 :
43 : namespace libMesh
44 : {
45 :
46 : //-----------------------------------------------------------------
47 : // JumpErrorEstimator implementations
48 0 : void JumpErrorEstimator::init_context (FEMContext &)
49 : {
50 : // Derived classes are *supposed* to rederive this
51 : libmesh_deprecated();
52 0 : }
53 :
54 :
55 :
56 25068 : void JumpErrorEstimator::estimate_error (const System & system,
57 : ErrorVector & error_per_cell,
58 : const NumericVector<Number> * solution_vector,
59 : bool estimate_parent_error)
60 : {
61 1616 : LOG_SCOPE("estimate_error()", "JumpErrorEstimator");
62 :
63 : /**
64 : * Conventions for assigning the direction of the normal:
65 : *
66 : * - e & f are global element ids
67 : *
68 : * Case (1.) Elements are at the same level, e<f
69 : * Compute the flux jump on the face and
70 : * add it as a contribution to error_per_cell[e]
71 : * and error_per_cell[f]
72 : *
73 : * ----------------------
74 : * | | |
75 : * | | f |
76 : * | | |
77 : * | e |---> n |
78 : * | | |
79 : * | | |
80 : * ----------------------
81 : *
82 : *
83 : * Case (2.) The neighbor is at a higher level.
84 : * Compute the flux jump on e's face and
85 : * add it as a contribution to error_per_cell[e]
86 : * and error_per_cell[f]
87 : *
88 : * ----------------------
89 : * | | | |
90 : * | | e |---> n |
91 : * | | | |
92 : * |-----------| f |
93 : * | | | |
94 : * | | | |
95 : * | | | |
96 : * ----------------------
97 : */
98 :
99 : // This parameter is not used when !LIBMESH_ENABLE_AMR.
100 808 : libmesh_ignore(estimate_parent_error);
101 :
102 : // The current mesh
103 1616 : const MeshBase & mesh = system.get_mesh();
104 :
105 : // The number of variables in the system
106 25068 : const unsigned int n_vars = system.n_vars();
107 :
108 : // The DofMap for this system
109 : #ifdef LIBMESH_ENABLE_AMR
110 808 : const DofMap & dof_map = system.get_dof_map();
111 : #endif
112 :
113 : // Resize the error_per_cell vector according to the
114 : // maximum element ID because we will be indexing it with IDs.
115 : // Initialize to 0.
116 25068 : error_per_cell.resize (mesh.max_elem_id());
117 808 : std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
118 :
119 : // Declare a vector of floats which is as long as
120 : // error_per_cell above, and fill with zeros. This vector will be
121 : // used to keep track of the number of edges (faces) on each active
122 : // element which are either:
123 : // 1) an internal edge
124 : // 2) an edge on a Neumann boundary for which a boundary condition
125 : // function has been specified.
126 : // The error estimator can be scaled by the number of flux edges (faces)
127 : // which the element actually has to obtain a more uniform measure
128 : // of the error. Use floats instead of ints since in case 2 (above)
129 : // f gets 1/2 of a flux face contribution from each of his
130 : // neighbors
131 1616 : std::vector<float> n_flux_faces;
132 25068 : if (scale_by_n_flux_faces)
133 0 : n_flux_faces.resize(error_per_cell.size(), 0);
134 :
135 : // Prepare current_local_solution to localize a non-standard
136 : // solution vector if necessary
137 25068 : if (solution_vector && solution_vector != system.solution.get())
138 : {
139 0 : NumericVector<Number> * newsol =
140 : const_cast<NumericVector<Number> *>(solution_vector);
141 0 : System & sys = const_cast<System &>(system);
142 0 : newsol->swap(*sys.solution);
143 0 : sys.update();
144 : }
145 :
146 : // We don't use full elem_jacobian or subjacobians here.
147 : fine_context = std::make_unique<FEMContext>
148 48520 : (system, nullptr, /* allocate_local_matrices = */ false);
149 : coarse_context = std::make_unique<FEMContext>
150 48520 : (system, nullptr, /* allocate_local_matrices = */ false);
151 :
152 : // Don't overintegrate - we're evaluating differences of FE values,
153 : // not products of them.
154 25068 : if (this->use_unweighted_quadrature_rules)
155 24125 : fine_context->use_unweighted_quadrature_rules(system.extra_quadrature_order);
156 :
157 : // Loop over all the variables we've been requested to find jumps in, to
158 : // pre-request
159 50136 : for (var=0; var<n_vars; var++)
160 : {
161 : // Skip variables which aren't part of our norm,
162 : // as well as SCALAR variables, which have no jumps
163 25876 : if (error_norm.weight(var) == 0.0 ||
164 25068 : system.variable_type(var).family == SCALAR)
165 0 : continue;
166 :
167 : // FIXME: Need to generalize this to vector-valued elements. [PB]
168 808 : FEBase * side_fe = nullptr;
169 :
170 : const std::set<unsigned char> & elem_dims =
171 808 : fine_context->elem_dimensions();
172 :
173 50278 : for (const auto & dim : elem_dims)
174 : {
175 25210 : fine_context->get_side_fe( var, side_fe, dim );
176 :
177 1978 : side_fe->get_xyz();
178 : }
179 : }
180 :
181 25876 : this->init_context(*fine_context);
182 25876 : this->init_context(*coarse_context);
183 :
184 : // If we're integrating jumps across mesh slits, then we'll need a
185 : // point locator to find slits, and we'll need to integrate point by
186 : // point on sides.
187 24260 : std::unique_ptr<PointLocatorBase> point_locator;
188 24260 : std::unique_ptr<const Elem> side_ptr;
189 :
190 25068 : if (integrate_slits)
191 138 : point_locator = mesh.sub_point_locator();
192 :
193 : // Iterate over all the active elements in the mesh
194 : // that live on this processor.
195 2208220 : for (const auto & e : mesh.active_local_element_ptr_range())
196 : {
197 1209949 : const dof_id_type e_id = e->id();
198 :
199 : #ifdef LIBMESH_ENABLE_AMR
200 :
201 376694 : if (e->infinite())
202 : {
203 : libmesh_warning("Warning: Jumps on the border of infinite elements are ignored."
204 : << std::endl);
205 1280 : continue;
206 : }
207 :
208 : // See if the parent of element e has been examined yet;
209 : // if not, we may want to compute the estimator on it
210 374774 : const Elem * parent = e->parent();
211 :
212 : // We only can compute and only need to compute on
213 : // parents with all active children
214 129863 : bool compute_on_parent = true;
215 1208029 : if (!parent || !estimate_parent_error)
216 129863 : compute_on_parent = false;
217 : else
218 0 : for (auto & child : parent->child_ref_range())
219 0 : if (!child.active())
220 0 : compute_on_parent = false;
221 :
222 129863 : if (compute_on_parent &&
223 0 : !error_per_cell[parent->id()])
224 : {
225 : // Compute a projection onto the parent
226 0 : DenseVector<Number> Uparent;
227 : FEBase::coarsened_dof_values
228 0 : (*(system.solution), dof_map, parent, Uparent, false);
229 :
230 : // Loop over the neighbors of the parent
231 0 : for (auto n_p : parent->side_index_range())
232 : {
233 0 : if (parent->neighbor_ptr(n_p) != nullptr) // parent has a neighbor here
234 : {
235 : // Find the active neighbors in this direction
236 0 : std::vector<const Elem *> active_neighbors;
237 : parent->neighbor_ptr(n_p)->
238 0 : active_family_tree_by_neighbor(active_neighbors,
239 : parent);
240 : // Compute the flux to each active neighbor
241 0 : for (std::size_t a=0,
242 0 : n_active_neighbors = active_neighbors.size();
243 0 : a != n_active_neighbors; ++a)
244 : {
245 0 : const Elem * f = active_neighbors[a];
246 :
247 0 : if (f ->infinite()) // don't take infinite elements into account
248 0 : continue;
249 :
250 : // FIXME - what about when f->level <
251 : // parent->level()??
252 0 : if (f->level() >= parent->level())
253 : {
254 0 : fine_context->pre_fe_reinit(system, f);
255 0 : coarse_context->pre_fe_reinit(system, parent);
256 0 : libmesh_assert_equal_to
257 : (coarse_context->get_elem_solution().size(),
258 : Uparent.size());
259 0 : coarse_context->get_elem_solution() = Uparent;
260 :
261 0 : this->reinit_sides();
262 :
263 : // Loop over all significant variables in the system
264 0 : for (var=0; var<n_vars; var++)
265 0 : if (error_norm.weight(var) != 0.0 &&
266 0 : system.variable_type(var).family != SCALAR)
267 : {
268 0 : this->internal_side_integration();
269 :
270 0 : error_per_cell[fine_context->get_elem().id()] +=
271 0 : static_cast<ErrorVectorReal>(fine_error);
272 0 : error_per_cell[coarse_context->get_elem().id()] +=
273 0 : static_cast<ErrorVectorReal>(coarse_error);
274 : }
275 :
276 : // Keep track of the number of internal flux
277 : // sides found on each element
278 0 : if (scale_by_n_flux_faces)
279 : {
280 0 : n_flux_faces[fine_context->get_elem().id()]++;
281 0 : n_flux_faces[coarse_context->get_elem().id()] +=
282 0 : this->coarse_n_flux_faces_increment();
283 : }
284 : }
285 : }
286 : }
287 0 : else if (integrate_boundary_sides)
288 : {
289 0 : fine_context->pre_fe_reinit(system, parent);
290 0 : libmesh_assert_equal_to
291 : (fine_context->get_elem_solution().size(),
292 : Uparent.size());
293 0 : fine_context->get_elem_solution() = Uparent;
294 0 : fine_context->side = cast_int<unsigned char>(n_p);
295 0 : fine_context->side_fe_reinit();
296 :
297 : // If we find a boundary flux for any variable,
298 : // let's just count it as a flux face for all
299 : // variables. Otherwise we'd need to keep track of
300 : // a separate n_flux_faces and error_per_cell for
301 : // every single var.
302 0 : bool found_boundary_flux = false;
303 :
304 0 : for (var=0; var<n_vars; var++)
305 0 : if (error_norm.weight(var) != 0.0 &&
306 0 : system.variable_type(var).family != SCALAR)
307 : {
308 0 : if (this->boundary_side_integration())
309 : {
310 0 : error_per_cell[fine_context->get_elem().id()] +=
311 0 : static_cast<ErrorVectorReal>(fine_error);
312 0 : found_boundary_flux = true;
313 : }
314 : }
315 :
316 0 : if (scale_by_n_flux_faces && found_boundary_flux)
317 0 : n_flux_faces[fine_context->get_elem().id()]++;
318 : }
319 : }
320 : }
321 : #endif // #ifdef LIBMESH_ENABLE_AMR
322 :
323 : // If we do any more flux integration, e will be the fine element
324 1208029 : fine_context->pre_fe_reinit(system, e);
325 :
326 : // Loop over the neighbors of element e
327 5602759 : for (auto n_e : e->side_index_range())
328 : {
329 4749043 : if ((e->neighbor_ptr(n_e) != nullptr) ||
330 112962 : integrate_boundary_sides)
331 : {
332 4151905 : fine_context->side = cast_int<unsigned char>(n_e);
333 4151905 : fine_context->side_fe_reinit();
334 : }
335 :
336 : // e is not on the boundary (infinite elements are treated as boundary)
337 4264867 : if (e->neighbor_ptr(n_e) != nullptr
338 4264867 : && !e->neighbor_ptr(n_e) ->infinite())
339 : {
340 :
341 1331074 : const Elem * f = e->neighbor_ptr(n_e);
342 921808 : const dof_id_type f_id = f->id();
343 :
344 : // Compute flux jumps if we are in case 1 or case 2.
345 4460268 : if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
346 2952246 : || (f->level() < e->level()))
347 : {
348 : // f is now the coarse element
349 2152201 : coarse_context->pre_fe_reinit(system, f);
350 :
351 2152201 : this->reinit_sides();
352 :
353 : // Loop over all significant variables in the system
354 4304402 : for (var=0; var<n_vars; var++)
355 2391149 : if (error_norm.weight(var) != 0.0 &&
356 2152201 : system.variable_type(var).family != SCALAR)
357 : {
358 2152201 : this->internal_side_integration();
359 :
360 2152201 : error_per_cell[fine_context->get_elem().id()] +=
361 2152201 : static_cast<ErrorVectorReal>(fine_error);
362 2152201 : error_per_cell[coarse_context->get_elem().id()] +=
363 2152201 : static_cast<ErrorVectorReal>(coarse_error);
364 : }
365 :
366 : // Keep track of the number of internal flux
367 : // sides found on each element
368 2152201 : if (scale_by_n_flux_faces)
369 : {
370 0 : n_flux_faces[fine_context->get_elem().id()]++;
371 0 : n_flux_faces[coarse_context->get_elem().id()] +=
372 0 : this->coarse_n_flux_faces_increment();
373 : }
374 : } // end if (case1 || case2)
375 : } // if (e->neighbor(n_e) != nullptr)
376 :
377 : // e might not have a neighbor_ptr, but might still have
378 : // another element sharing its side. This can happen in a
379 : // mesh where solution continuity is maintained via nodal
380 : // constraint rows.
381 115332 : else if (integrate_slits)
382 : {
383 832 : side_ptr = e->build_side_ptr(n_e);
384 64 : std::set<const Elem *> candidate_elements;
385 768 : (*point_locator)(side_ptr->vertex_average(), candidate_elements);
386 :
387 : // We should have at least found ourselves...
388 64 : libmesh_assert(candidate_elements.count(e));
389 :
390 : // If we only found ourselves, this probably isn't a
391 : // slit; we don't yet support meshes so non-conforming
392 : // as to have overlap of part of an element side without
393 : // overlap of its center.
394 768 : if (candidate_elements.size() < 2)
395 48 : continue;
396 :
397 192 : FEType hardest_fe_type = fine_context->find_hardest_fe_type();
398 :
399 192 : auto dim = e->dim();
400 :
401 : auto side_qrule =
402 : hardest_fe_type.unweighted_quadrature_rule
403 208 : (dim-1, system.extra_quadrature_order);
404 208 : auto side_fe = FEAbstract::build(dim, hardest_fe_type);
405 208 : side_fe->attach_quadrature_rule(side_qrule.get());
406 48 : const std::vector<Point> & qface_point = side_fe->get_xyz();
407 192 : side_fe->reinit(e, n_e);
408 :
409 384 : for (auto qp : make_range(side_qrule->n_points()))
410 : {
411 192 : const Point p = qface_point[qp];
412 208 : const std::vector<Point> qp_pointvec(1, p);
413 32 : std::set<const Elem *> side_elements;
414 192 : (*point_locator)(side_ptr->vertex_average(), side_elements);
415 :
416 : // If we have multiple neighbors meeting here we'll just
417 : // take weighted jumps from all of them.
418 : //
419 : // We'll also do integrations from both sides of slits,
420 : // rather than try to figure out a disambiguation rule
421 : // that makes sense for non-conforming slits in general.
422 : // This means we want an extra factor of 0.5 on the
423 : // integrals to compensate for doubling them.
424 192 : const std::size_t n_neighbors = side_elements.size() - 1;
425 192 : const Real neighbor_frac = Real(1)/n_neighbors;
426 :
427 : const std::vector<Real>
428 208 : qp_weightvec(1, neighbor_frac * side_qrule->w(qp));
429 :
430 576 : for (const Elem * f : side_elements)
431 : {
432 384 : if (f == e)
433 192 : continue;
434 :
435 192 : coarse_context->pre_fe_reinit(system, f);
436 192 : fine_context->pre_fe_reinit(system, e);
437 32 : std::vector<Point> qp_coarse, qp_fine;
438 384 : for (unsigned int v=0; v<n_vars; v++)
439 208 : if (error_norm.weight(v) != 0.0 &&
440 192 : fine_context->get_system().variable_type(v).family != SCALAR)
441 : {
442 16 : FEBase * coarse_fe = coarse_context->get_side_fe(v, dim);
443 192 : if (qp_coarse.empty())
444 192 : FEMap::inverse_map (dim, f, qp_pointvec, qp_coarse);
445 192 : coarse_fe->reinit(f, &qp_coarse, &qp_weightvec);
446 16 : FEBase * fine_fe = fine_context->get_side_fe(v, dim);
447 192 : if (qp_fine.empty())
448 192 : FEMap::inverse_map (dim, e, qp_pointvec, qp_fine);
449 192 : fine_fe->reinit(e, &qp_fine, &qp_weightvec);
450 : }
451 :
452 : // Loop over all significant variables in the system
453 384 : for (var=0; var<n_vars; var++)
454 208 : if (error_norm.weight(var) != 0.0 &&
455 192 : system.variable_type(var).family != SCALAR)
456 : {
457 192 : this->internal_side_integration();
458 :
459 192 : error_per_cell[fine_context->get_elem().id()] +=
460 192 : static_cast<ErrorVectorReal>(fine_error);
461 192 : error_per_cell[coarse_context->get_elem().id()] +=
462 192 : static_cast<ErrorVectorReal>(coarse_error);
463 : }
464 : }
465 : }
466 160 : }
467 :
468 : // Otherwise, e is on the boundary. If it happens to
469 : // be on a Dirichlet boundary, we need not do anything.
470 : // On the other hand, if e is on a Neumann (flux) boundary
471 : // with grad(u).n = g, we need to compute the additional residual
472 : // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
473 : // We can only do this with some knowledge of the boundary
474 : // conditions, i.e. the user must have attached an appropriate
475 : // BC function.
476 114564 : else if (integrate_boundary_sides)
477 : {
478 0 : if (integrate_slits)
479 0 : libmesh_not_implemented();
480 :
481 0 : bool found_boundary_flux = false;
482 :
483 0 : for (var=0; var<n_vars; var++)
484 0 : if (error_norm.weight(var) != 0.0 &&
485 0 : system.variable_type(var).family != SCALAR)
486 0 : if (this->boundary_side_integration())
487 : {
488 0 : error_per_cell[fine_context->get_elem().id()] +=
489 0 : static_cast<ErrorVectorReal>(fine_error);
490 0 : found_boundary_flux = true;
491 : }
492 :
493 0 : if (scale_by_n_flux_faces && found_boundary_flux)
494 0 : n_flux_faces[fine_context->get_elem().id()]++;
495 : } // end if (e->neighbor_ptr(n_e) == nullptr)
496 : } // end loop over neighbors
497 23452 : } // End loop over active local elements
498 :
499 :
500 : // Each processor has now computed the error contributions
501 : // for its local elements. We need to sum the vector
502 : // and then take the square-root of each component. Note
503 : // that we only need to sum if we are running on multiple
504 : // processors, and we only need to take the square-root
505 : // if the value is nonzero. There will in general be many
506 : // zeros for the inactive elements.
507 :
508 : // First sum the vector of estimated error values
509 25068 : this->reduce_error(error_per_cell, system.comm());
510 :
511 : // Compute the square-root of each component.
512 9062775 : for (auto i : index_range(error_per_cell))
513 9384173 : if (error_per_cell[i] != 0.)
514 7046792 : error_per_cell[i] = std::sqrt(error_per_cell[i]);
515 :
516 :
517 25068 : if (this->scale_by_n_flux_faces)
518 : {
519 : // Sum the vector of flux face counts
520 0 : this->reduce_error(n_flux_faces, system.comm());
521 :
522 : // Sanity check: Make sure the number of flux faces is
523 : // always an integer value
524 : #ifdef DEBUG
525 0 : for (const auto & val : n_flux_faces)
526 0 : libmesh_assert_equal_to (val, static_cast<float>(static_cast<unsigned int>(val)));
527 : #endif
528 :
529 : // Scale the error by the number of flux faces for each element
530 0 : for (auto i : index_range(n_flux_faces))
531 : {
532 0 : if (n_flux_faces[i] == 0.0) // inactive or non-local element
533 0 : continue;
534 :
535 0 : error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
536 : }
537 : }
538 :
539 : // If we used a non-standard solution before, now is the time to fix
540 : // the current_local_solution
541 25068 : if (solution_vector && solution_vector != system.solution.get())
542 : {
543 0 : NumericVector<Number> * newsol =
544 : const_cast<NumericVector<Number> *>(solution_vector);
545 0 : System & sys = const_cast<System &>(system);
546 0 : newsol->swap(*sys.solution);
547 0 : sys.update();
548 : }
549 48520 : }
550 :
551 :
552 :
553 : void
554 2152201 : JumpErrorEstimator::reinit_sides ()
555 : {
556 2152201 : fine_context->side_fe_reinit();
557 :
558 2152201 : unsigned short dim = fine_context->get_elem().dim();
559 238948 : libmesh_assert_equal_to(dim, coarse_context->get_elem().dim());
560 :
561 238948 : FEBase * fe_fine = nullptr;
562 238948 : fine_context->get_side_fe( 0, fe_fine, dim );
563 :
564 : // Get the physical locations of the fine element quadrature points
565 2391149 : std::vector<Point> qface_point = fe_fine->get_xyz();
566 :
567 : // Find the master quadrature point locations on the coarse element
568 238948 : FEBase * fe_coarse = nullptr;
569 238948 : coarse_context->get_side_fe( 0, fe_coarse, dim );
570 :
571 477896 : std::vector<Point> qp_coarse;
572 :
573 2152201 : FEMap::inverse_map (coarse_context->get_elem().dim(),
574 477896 : &coarse_context->get_elem(), qface_point,
575 : qp_coarse);
576 :
577 : // The number of variables in the system
578 238948 : const unsigned int n_vars = fine_context->n_vars();
579 :
580 : // Calculate all coarse element shape functions at those locations
581 4304402 : for (unsigned int v=0; v<n_vars; v++)
582 2391149 : if (error_norm.weight(v) != 0.0 &&
583 2152201 : fine_context->get_system().variable_type(v).family != SCALAR)
584 : {
585 238948 : coarse_context->get_side_fe( v, fe_coarse, dim );
586 2152201 : fe_coarse->reinit (&coarse_context->get_elem(), &qp_coarse);
587 : }
588 2152201 : }
589 :
590 :
591 :
592 0 : float JumpErrorEstimator::coarse_n_flux_faces_increment ()
593 : {
594 : // Keep track of the number of internal flux sides found on each
595 : // element
596 0 : unsigned short dim = coarse_context->get_elem().dim();
597 :
598 : const unsigned int divisor =
599 0 : 1 << (dim-1)*(fine_context->get_elem().level() -
600 0 : coarse_context->get_elem().level());
601 :
602 : // With a difference of n levels between fine and coarse elements,
603 : // we compute a fractional flux face for the coarse element by adding:
604 : // 1/2^n in 2D
605 : // 1/4^n in 3D
606 : // each time. This code will get hit 2^n times in 2D and 4^n
607 : // times in 3D so that the final flux face count for the coarse
608 : // element will be an integer value.
609 :
610 0 : return 1.0f / static_cast<float>(divisor);
611 : }
612 :
613 : } // namespace libMesh
|
