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