libMesh/libmesh: coverage diff

Base 82cc40 Head #4232 290bfc
Total Total +/- New
Rate 64.57% 64.72% +0.14% 96.97%
Hits 75998 76189 +191 32
Misses 41693 41534 -159 1
Filename Stmts Miss Cover
include/systems/variational_smoother_constraint.h 0 -1 +3.85%
src/geom/cell_hex20.C 0 +2 -0.71%
src/geom/cell_hex27.C 0 +2 -0.47%
src/geom/cell_pyramid13.C 0 -77 +36.32%
src/geom/cell_pyramid14.C 0 -80 +33.90%
src/mesh/mesh_refinement.C 0 -2 +0.29%
src/mesh/mesh_refinement_smoothing.C 0 -1 +0.39%
src/systems/variational_smoother_constraint.C 0 -3 +0.80%
src/systems/variational_smoother_system.C +32 +1 -0.12%
TOTAL +32 -159 +0.14%
code
coverage unchanged
code
coverage increased
code
coverage decreased
+
line added or modified

include/systems/variational_smoother_constraint.h

81  
82  
83  
84  
85  
86  
87  

   * @param p The point in question
   * @return bool indicating whether p lies on this point.
   */
  bool contains_point(const PointConstraint &p) const { return *this == p; }

  /**
   * Computes the intersection of this point with another constraint.

src/geom/cell_hex20.C

140  
141  
142  
143  
144  
145  
146  
147  
148  
149  
150  
151  
152  
153  

      !v.relative_fuzzy_equals(this->point(13) - this->point(1), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(14) - this->point(2), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(15) - this->point(3), affine_tol))
    return false;
  // Make sure y-edges are straight
  v = (this->point(3) - this->point(0))/2;
  if (!v.relative_fuzzy_equals(this->point(11) - this->point(0), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(9) - this->point(1), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(17) - this->point(5), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(19) - this->point(4), affine_tol))
    return false;
  // If all the above checks out, the map is affine
  return true;
}

src/geom/cell_hex27.C

154  
155  
156  
157  
158  
159  
160  
161  
162  
163  
164  
165  
166  
167  

      !v.relative_fuzzy_equals(this->point(15) - this->point(3), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(22) - this->point(9), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(24) - this->point(11), affine_tol))
    return false;
  // Make sure y-edges are straight
  v = (this->point(3) - this->point(0))/2;
  if (!v.relative_fuzzy_equals(this->point(11) - this->point(0), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(9) - this->point(1), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(17) - this->point(5), affine_tol) ||
      !v.relative_fuzzy_equals(this->point(19) - this->point(4), affine_tol))
    return false;
  // If all the above checks out, the map is affine
  return true;
}

src/geom/cell_pyramid13.C

349  
350  
351  
352  
353  
354  
355  
356  
357  
358  
359  
360  
361  
362  
363  
364  
365  
366  
367  
368  
369  
370  
371  
372  
373  
374  
375  
376  
377  
378  
379  
380  
381  
382  
383  
384  
385  
386  
387  
388  
389  
390  
391  
392  
393  
394  
395  
396  
397  
398  
399  
400  
401  
402  
403  
404  
405  
406  
407  
408  
409  
410  
411  
412  
413  
414  
415  
416  
417  
418  
419  
420  
421  
422  
423  
424  
425  
426  
427  
428  
429  
430  




Real Pyramid13::volume () const
{
  // This specialization is good for Lagrange mappings only
  if (this->mapping_type() != LAGRANGE_MAP)
    return this->Elem::volume();

  // Make copies of our points.  It makes the subsequent calculations a bit
  // shorter and avoids dereferencing the same pointer multiple times.
  Point
    x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),   x4 = point(4),   x5 = point(5),   x6 = point(6),
    x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11), x12 = point(12);

  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
  // These are copied directly from the output of a Python script.
  Point dx_dxi[14] =
    {
      x6/2 - x8/2,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
      -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
      x0/4 - x1/4 + x2/4 - x3/4,
      -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
      x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
      -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
      x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
      x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
      -x0 - x1 - x2 - x3 + 2*x5 + 2*x7,
      x0/2 + x1/2 + x2/2 + x3/2 - x5 - x7,
      -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
      x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7
    };

  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
  // These are copied directly from the output of a Python script.
  Point dx_deta[14] =
    {
      -x5/2 + x7/2,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
      -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
      x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
      -x0 - x1 - x2 - x3 + 2*x6 + 2*x8,
      x0/2 + x1/2 + x2/2 + x3/2 - x6 - x8,
      x0/4 - x1/4 + x2/4 - x3/4,
      -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
      x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
      -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
      x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
      -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
      x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
    };

  // dx/dxi and dx/deta have 14 components while dx/dzeta has 19.
  // These are copied directly from the output of a Python script.
  Point dx_dzeta[19] =
    {
      x0/4 + x1/4 + x10 + x11 + x12 + x2/4 + x3/4 - x4 - x5 - x6 - x7 - x8 + x9,
      -3*x0/4 - 3*x1/4 - 5*x10 - 5*x11 - 5*x12 - 3*x2/4 - 3*x3/4 + 7*x4 + 4*x5 + 4*x6 + 4*x7 + 4*x8 - 5*x9,
      3*x0/4 + 3*x1/4 + 9*x10 + 9*x11 + 9*x12 + 3*x2/4 + 3*x3/4 - 15*x4 - 6*x5 - 6*x6 - 6*x7 - 6*x8 + 9*x9,
      -x0/4 - x1/4 - 7*x10 - 7*x11 - 7*x12 - x2/4 - x3/4 + 13*x4 + 4*x5 + 4*x6 + 4*x7 + 4*x8 - 7*x9,
      2*x10 + 2*x11 + 2*x12 - 4*x4 - x5 - x6 - x7 - x8 + 2*x9,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
      -3*x0/4 - 3*x1/4 + 3*x10 - 3*x11 - 3*x12 + 3*x2/4 + 3*x3/4 - 3*x5/2 + 3*x7/2 + 3*x9,
      3*x0/4 + 3*x1/4 - 3*x10 + 3*x11 + 3*x12 - 3*x2/4 - 3*x3/4 + 3*x5/2 - 3*x7/2 - 3*x9,
      -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
      -3*x0/4 + 3*x1/4 - 3*x10 - 3*x11 + 3*x12 + 3*x2/4 - 3*x3/4 + 3*x6/2 - 3*x8/2 + 3*x9,
      3*x0/4 - 3*x1/4 + 3*x10 + 3*x11 - 3*x12 - 3*x2/4 + 3*x3/4 - 3*x6/2 + 3*x8/2 - 3*x9,
      -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
      -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
      x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
      -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
      x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
      -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
      x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2
    };

  // The (xi, eta, zeta) exponents for each of the dx_dxi terms
  static const int dx_dxi_exponents[14][3] =

489  
490  
491  
492  
493  
494  
495  

    };

  // Number of points in the quadrature rule
  const int N = 27;

  // Parameters of the quadrature rule
  static const Real

613  
614  
615  
616  
617  
618  
619  
620  
621  
622  
623  
624  
625  
626  
627  
628  
629  
630  
631  
632  
633  
634  
635  
636  
637  
638  
639  
640  
641  
642  
643  
644  
645  
646  
647  
648  
649  
650  
651  
652  
653  
654  
655  
656  
657  
658  
659  
660  
661  

                            w1, w2, w3, w4, w5, w6, // 18-23
                            w1, w2, w3};            // 24-26

  Real vol = 0.;
  for (int q=0; q<N; ++q)
    {
      // Compute denominators for the current q.
      Real
        den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
        den3 = den2*(1. - zeta[q][1]);

      // Compute dx/dxi and dx/deta at the current q.
      Point dx_dxi_q, dx_deta_q;
      for (int c=0; c<14; ++c)
        {
          dx_dxi_q +=
            xi[q][dx_dxi_exponents[c][0]]*
            eta[q][dx_dxi_exponents[c][1]]*
            zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];

          dx_deta_q +=
            xi[q][dx_deta_exponents[c][0]]*
            eta[q][dx_deta_exponents[c][1]]*
            zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
        }

      // Compute dx/dzeta at the current q.
      Point dx_dzeta_q;
      for (int c=0; c<19; ++c)
        {
          dx_dzeta_q +=
            xi[q][dx_dzeta_exponents[c][0]]*
            eta[q][dx_dzeta_exponents[c][1]]*
            zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
        }

      // Scale everything appropriately
      dx_dxi_q /= den2;
      dx_deta_q /= den2;
      dx_dzeta_q /= den3;

      // Compute scalar triple product, multiply by weight, and accumulate volume.
      vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
    }

  return vol;
}

src/geom/cell_pyramid14.C

387  
388  
389  
390  
391  
392  
393  
394  
395  
396  
397  
398  
399  
400  
401  
402  
403  
404  
405  
406  
407  
408  
409  
410  
411  
412  
413  
414  
415  
416  
417  
418  
419  
420  
421  
422  
423  
424  
425  
426  
427  
428  
429  
430  
431  
432  
433  
434  
435  
436  
437  
438  
439  
440  
441  
442  
443  
444  
445  
446  
447  
448  
449  
450  
451  
452  
453  
454  
455  
456  
457  
458  
459  
460  
461  
462  
463  
464  
465  
466  
467  
468  
469  
470  
471  




Real Pyramid14::volume () const
{
  // This specialization is good for Lagrange mappings only
  if (this->mapping_type() != LAGRANGE_MAP)
    return this->Elem::volume();

  // Make copies of our points.  It makes the subsequent calculations a bit
  // shorter and avoids dereferencing the same pointer multiple times.
  Point
    x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),   x4 = point(4),   x5 = point(5),   x6 = point(6),
    x7 = point(7), x8 = point(8), x9 = point(9), x10 = point(10), x11 = point(11), x12 = point(12), x13 = point(13);

  // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
  // These are copied directly from the output of a Python script.
  Point dx_dxi[15] =
    {
      x6/2 - x8/2,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - 3*x6/2 + 3*x8/2 - x9,
      -x0/2 + x1/2 - 2*x10 - 2*x11 + 2*x12 + x2/2 - x3/2 + 3*x6/2 - 3*x8/2 + 2*x9,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
      x0/4 - x1/4 + x2/4 - x3/4,
      -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
      x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
      -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
      x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
      -2*x13 + x6 + x8,
      4*x13 - 2*x6 - 2*x8,
      -2*x13 + x6 + x8,
      -x0/2 - x1/2 + x2/2 + x3/2 + x5 - x7,
      x0/2 + x1/2 - x2/2 - x3/2 - x5 + x7,
      x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
    };

  // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
  // These are copied directly from the output of a Python script.
  Point dx_deta[15] =
    {
      -x5/2 + x7/2,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + 3*x5/2 - 3*x7/2 - x9,
      -x0/2 - x1/2 + 2*x10 - 2*x11 - 2*x12 + x2/2 + x3/2 - 3*x5/2 + 3*x7/2 + 2*x9,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
      -2*x13 + x5 + x7,
      4*x13 - 2*x5 - 2*x7,
      -2*x13 + x5 + x7,
      x0/4 - x1/4 + x2/4 - x3/4,
      -3*x0/4 + 3*x1/4 - x10 + x11 - x12 - 3*x2/4 + 3*x3/4 + x9,
      x0/2 - x1/2 + x10 - x11 + x12 + x2/2 - x3/2 - x9,
      -x0/2 + x1/2 + x2/2 - x3/2 - x6 + x8,
      x0/2 - x1/2 - x2/2 + x3/2 + x6 - x8,
      -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
      x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2,
      x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
    };

  // dx/dxi and dx/deta have 15 components while dx/dzeta has 20.
  // These are copied directly from the output of a Python script.
  Point dx_dzeta[20] =
    {
      -x0/4 - x1/4 + x10 + x11 + x12 - 2*x13 - x2/4 - x3/4 - x4 + x9,
      5*x0/4 + 5*x1/4 - 5*x10 - 5*x11 - 5*x12 + 8*x13 + 5*x2/4 + 5*x3/4 + 7*x4 - 5*x9,
      -9*x0/4 - 9*x1/4 + 9*x10 + 9*x11 + 9*x12 - 12*x13 - 9*x2/4 - 9*x3/4 - 15*x4 + 9*x9,
      7*x0/4 + 7*x1/4 - 7*x10 - 7*x11 - 7*x12 + 8*x13 + 7*x2/4 + 7*x3/4 + 13*x4 - 7*x9,
      -x0/2 - x1/2 + 2*x10 + 2*x11 + 2*x12 - 2*x13 - x2/2 - x3/2 - 4*x4 + 2*x9,
      x0/4 + x1/4 - x10 + x11 + x12 - x2/4 - x3/4 + x5/2 - x7/2 - x9,
      -3*x0/4 - 3*x1/4 + 3*x10 - 3*x11 - 3*x12 + 3*x2/4 + 3*x3/4 - 3*x5/2 + 3*x7/2 + 3*x9,
      3*x0/4 + 3*x1/4 - 3*x10 + 3*x11 + 3*x12 - 3*x2/4 - 3*x3/4 + 3*x5/2 - 3*x7/2 - 3*x9,
      -x0/4 - x1/4 + x10 - x11 - x12 + x2/4 + x3/4 - x5/2 + x7/2 + x9,
      x0/4 - x1/4 + x10 + x11 - x12 - x2/4 + x3/4 - x6/2 + x8/2 - x9,
      -3*x0/4 + 3*x1/4 - 3*x10 - 3*x11 + 3*x12 + 3*x2/4 - 3*x3/4 + 3*x6/2 - 3*x8/2 + 3*x9,
      3*x0/4 - 3*x1/4 + 3*x10 + 3*x11 - 3*x12 - 3*x2/4 + 3*x3/4 - 3*x6/2 + 3*x8/2 - 3*x9,
      -x0/4 + x1/4 - x10 - x11 + x12 + x2/4 - x3/4 + x6/2 - x8/2 + x9,
      -x0/4 + x1/4 - x10 + x11 - x12 - x2/4 + x3/4 + x9,
      x0/4 - x1/4 + x10 - x11 + x12 + x2/4 - x3/4 - x9,
      -x0/4 + x1/4 + x2/4 - x3/4 - x6/2 + x8/2,
      x0/4 - x1/4 - x2/4 + x3/4 + x6/2 - x8/2,
      -x0/4 - x1/4 + x2/4 + x3/4 + x5/2 - x7/2,
      x0/4 + x1/4 - x2/4 - x3/4 - x5/2 + x7/2,
      x0/2 + x1/2 + 2*x13 + x2/2 + x3/2 - x5 - x6 - x7 - x8
    };

  // The (xi, eta, zeta) exponents for each of the dx_dxi terms
  static const int dx_dxi_exponents[15][3] =

533  
534  
535  
536  
537  
538  
539  

    };

  // Number of points in the quadrature rule
  const int N = 27;

  // Parameters of the quadrature rule
  static const Real

657  
658  
659  
660  
661  
662  
663  
664  
665  
666  
667  
668  
669  
670  
671  
672  
673  
674  
675  
676  
677  
678  
679  
680  
681  
682  
683  
684  
685  
686  
687  
688  
689  
690  
691  
692  
693  
694  
695  
696  
697  
698  
699  
700  
701  
702  
703  
704  
705  

                            w1, w2, w3, w4, w5, w6, // 18-23
                            w1, w2, w3};            // 24-26

  Real vol = 0.;
  for (int q=0; q<N; ++q)
    {
      // Compute denominators for the current q.
      Real
        den2 = (1. - zeta[q][1])*(1. - zeta[q][1]),
        den3 = den2*(1. - zeta[q][1]);

      // Compute dx/dxi and dx/deta at the current q.
      Point dx_dxi_q, dx_deta_q;
      for (int c=0; c<15; ++c)
        {
          dx_dxi_q +=
            xi[q][dx_dxi_exponents[c][0]]*
            eta[q][dx_dxi_exponents[c][1]]*
            zeta[q][dx_dxi_exponents[c][2]]*dx_dxi[c];

          dx_deta_q +=
            xi[q][dx_deta_exponents[c][0]]*
            eta[q][dx_deta_exponents[c][1]]*
            zeta[q][dx_deta_exponents[c][2]]*dx_deta[c];
        }

      // Compute dx/dzeta at the current q.
      Point dx_dzeta_q;
      for (int c=0; c<20; ++c)
        {
          dx_dzeta_q +=
            xi[q][dx_dzeta_exponents[c][0]]*
            eta[q][dx_dzeta_exponents[c][1]]*
            zeta[q][dx_dzeta_exponents[c][2]]*dx_dzeta[c];
        }

      // Scale everything appropriately
      dx_dxi_q /= den2;
      dx_deta_q /= den2;
      dx_dzeta_q /= den3;

      // Compute scalar triple product, multiply by weight, and accumulate volume.
      vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
    }

  return vol;
}

src/mesh/mesh_refinement.C

1272  
1273  
1274  
1275  
1276  
1277  
1278  
1279  

                              if (neighbor->p_level() < my_p_level &&
                                  neighbor->p_refinement_flag() != Elem::REFINE)
                                {
                                  neighbor->set_p_refinement_flag(Elem::REFINE);
                                  level_one_satisfied = false;
                                  compatible_with_coarsening = false;
                                }
                              if (neighbor->p_level() == my_p_level &&

src/mesh/mesh_refinement_smoothing.C

531  
532  
533  
534  
535  
536  
537  

          if (elem->p_refinement_flag() == Elem::COARSEN)
            elem->set_p_refinement_flag(Elem::DO_NOTHING);
          else
            elem->set_p_refinement_flag(Elem::REFINE);
          flags_changed = true;
        }
    }

src/systems/variational_smoother_constraint.C

57  
58  
59  
60  
61  
62  
63  
64  
65  

  return _point < other.point();
}

bool PointConstraint::operator==(const PointConstraint & other) const
{
  return _point.absolute_fuzzy_equals(other.point(), _tol);
}

ConstraintVariant PointConstraint::intersect(const ConstraintVariant & other) const

274  
275  
276  
277  
278  
279  
280  

        else if constexpr (std::is_same_v<T, LineConstraint>)
          {
            if (this->contains_line(o))
              return o;

            if (this->is_parallel(o))
              // Line is parallel and does not intersect the plane

src/systems/variational_smoother_system.C

928  
929  
930  
931 +
932  
933  
934  



  // Elems deriving from Pyramid
  else if (type_str.compare(0, 7, "PYRAMID") == 0)
    {

      // The target element is a pyramid with an square base and

941  
942  
943  
944 +
945  
946 +
947  
948 +
949  
950 +
951 +
952  
953  
954 +
955 +
956 +
957 +
958 +
959  
960 +
961  
962 +
963  
964  
965  
966 +
967 +
968 +
969 +
970  
971  
972 +
973 +
974 +
975 +
976  
977 +
978  
979  
980 +
981 +
982  
983 +
984  
985  
986 +
987 +
988 +
989 +
990  
991 +
992  
993  
994 +
995 +
996  
997  
998  

      // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
      // non-optimal reference element.

      const Real sqrt_2 = std::sqrt(Real(2));
      // Side length that preserves the volume of the reference element
      const auto side_length = std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
      // Pyramid height with the property that all faces are equilateral triangles
      const auto target_height = side_length / sqrt_2;

      const auto & s = side_length;
      const auto & h = target_height;

      //                                         x        y        z    node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,      0.), 0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,      0.), 1));
      owned_nodes.emplace_back(Node::build(Point(s,       s,       0.), 2));
      owned_nodes.emplace_back(Node::build(Point(0.,      s,       0.), 3));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h),  4));

      if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
        {
          const auto & on = owned_nodes;
          // Define the edge midpoint nodes of the pyramid

          // Base node to base node midpoint nodes
          owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
          owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
          owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
          owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));

          // Base node to apex node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));

          if (type == PYRAMID14 || type == PYRAMID18)
            {
              // Define the square face midpoint node of the pyramid
              owned_nodes.emplace_back(
                  Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));

              if (type == PYRAMID18)
                {
                  // Define the triangular face nodes
                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
                  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
                  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
                }
            }
        }

      else if (type != PYRAMID5)
        libmesh_error_msg("Unsupported pyramid element: " << type_str);

    } // if Pyramid


1005  
1006  
1007  
1008 +
1009  
1010  
1011  

  for (const auto & node_ptr : owned_nodes)
    target_elem->set_node(node_ptr->id(), node_ptr.get());

  libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));

  return std::make_pair(std::move(target_elem), std::move(owned_nodes));
}