mesh_extras.hpp

Go to the documentation of this file.

// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
 
#ifndef MFEM_MESH_EXTRAS
#define MFEM_MESH_EXTRAS
 
#include "mfem.hpp"
#include <sstream>
 
namespace mfem
{
 
namespace common
{
 
class ElementMeshStream : public std::stringstream
{
public:
   ElementMeshStream(Element::Type e);
};
 
/// Merges vertices which lie at the same location
void MergeMeshNodes(Mesh * mesh, int logging);
 
/// Convert a set of attribute numbers to a marker array
/** The marker array will be of size max_attr and it will contain only zeroes
    and ones. Ones indicate which attribute numbers are present in the attrs
    array. In the special case when attrs has an entry equal to -1 the marker
    array will contain all ones. */
inline
void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
{
   if (attrs.Find(-1) != -1) { (marker = Array<int>(max_attr)) = 1; }
   else { marker = AttributeSets::AttrToMarker(max_attr, attrs); }
}
 
/// Transform a mesh according to an arbitrary affine transformation
///    y = A x + b
/// Where A is a spaceDim x spaceDim matrix and b is a vector of size spaceDim.
/// If A is of size zero the transformation will be y = b.
/// If b is of size zero the transformation will be y = A x.
///
/// Note that no error checking related to the determinant of A is performed.
/// If A has a non-positive determinant it is likely to produce an invalid
/// transformed mesh.
class AffineTransformation : public VectorCoefficient
{
private:
   DenseMatrix A;
   Vector b;
   Vector x;
 
public:
   AffineTransformation(int dim_, const DenseMatrix &A_, const Vector & b_)
      : VectorCoefficient(dim_), A(A_), b(b_), x(dim_)
   {
      MFEM_VERIFY((A.Height() == dim_ && A.Width() == dim_) ||
                  (A.Height() == 0 && A.Width() == 0),
                  "Affine transformation given an invalid matrix");
      MFEM_VERIFY(b.Size() == dim_ || b.Size() == 0,
                  "Affine transformation given an invalid vector");
   }
 
   void Eval(Vector &V, ElementTransformation &T,
             const IntegrationPoint &ip) override;
 
   using VectorCoefficient::Eval;
};
 
/// Generalized Kershaw mesh transformation in 2D and 3D, see D. Kershaw,
/// "Differencing of the diffusion equation in Lagrangian hydrodynamic codes",
/// JCP, 39:375–395, 1981.
/** The input mesh should be Cartesian nx x ny x nz with nx divisible by 6 and
    ny, nz divisible by 2.
    The parameters @a epsy and @a epsz must be in (0, 1].
    Uniform mesh is recovered for epsy=epsz=1.
    The @a smooth parameter controls the transition between different layers. */
// Usage:
// common::KershawTransformation kershawT(pmesh->Dimension(), 0.3, 0.3, 2);
// pmesh->Transform(kershawT);
class KershawTransformation : public VectorCoefficient
{
private:
   int dim;
   real_t epsy, epsz;
   int smooth;
 
public:
   KershawTransformation(const int dim_, real_t epsy_ = 0.3,
                         real_t epsz_ = 0.3, int smooth_ = 1)
      : VectorCoefficient(dim_), dim(dim_), epsy(epsy_),
        epsz(epsz_), smooth(smooth_)
   {
      MFEM_VERIFY(dim > 1,"Kershaw transformation only works for 2D and 3D"
                  "meshes.");
      MFEM_VERIFY(smooth >= 1 && smooth <= 3,
                  "Kershaw parameter smooth must be in [1, 3]");
      MFEM_VERIFY(epsy > 0 && epsy <=1,
                  "Kershaw parameter epsy must be in (0, 1].");
      if (dim == 3)
      {
         MFEM_VERIFY(epsz > 0 && epsz <=1,
                     "Kershaw parameter epsz must be in (0, 1].");
      }
   }
 
   // 1D transformation at the right boundary.
   real_t right(const real_t eps, const real_t x)
   {
      return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1);
   }
 
   // 1D transformation at the left boundary
   real_t left(const real_t eps, const real_t x)
   {
      return 1-right(eps,1-x);
   }
 
   // Transition from a value of "a" for x=0, to a value of "b" for x=1.
   // Controlled through "smooth" parameter.
   real_t step(const real_t a, const real_t b, real_t x)
   {
      if (x <= 0) { return a; }
      if (x >= 1) { return b; }
      if (smooth == 1) { return a + (b-a) * (x); }
      else if (smooth == 2) { return a + (b-a) * (x*x*(3-2*x)); }
      else { return a + (b-a) * (x*x*x*(x*(6*x-15)+10)); }
   }
 
   void Eval(Vector &V, ElementTransformation &T,
             const IntegrationPoint &ip) override;
 
   using VectorCoefficient::Eval;
};
 
/// Transform a [0,1]^D mesh into a spiral. The parameters are:
/// @a turns - number of turns around the origin,
/// @a width - for D >= 2, the width of the spiral arm,
/// @ gap    - gap between adjacent spiral arms at the end of each turn,
/// @ height - for D = 3, the maximum height of the spiral.
// Usage:
// common::SpiralTransformation spiralT(spaceDim, 2.4, 0.1, 0.05, 1.0);
// pmesh->Transform(spiralT);
class SpiralTransformation : public VectorCoefficient
{
private:
   real_t dim, turns, width, gap, height;
 
public:
   SpiralTransformation(int dim_, real_t turns_ = 1.0, real_t width_ = 0.1,
                        real_t gap_ = 0.05, real_t height_ = 1.0)
      : VectorCoefficient(dim_), dim(dim_),
        turns(turns_), width(width_), gap(gap_), height(height_)
   {
      MFEM_VERIFY(turns > 0 && width > 0 && gap > 0 && height > 0,
                  "Spiral transformation requires positive parameters: turns, "
                  " width, gap, and height.");
   }
 
   void Eval(Vector &V, ElementTransformation &T,
             const IntegrationPoint &ip) override;
 
   using VectorCoefficient::Eval;
};
 
 
} // namespace common
 
} // namespace mfem
 
#endif