lor_mms.hpp

Go to the documentation of this file.

// Copyright (c) 2010-2024, 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_LOR_MMS_HPP
#define MFEM_LOR_MMS_HPP
 
namespace mfem
{
 
static constexpr real_t pi = M_PI, pi2 = M_PI*M_PI;
 
// Exact solution for definite Helmholtz problem with RHS corresponding to f
// defined below.
real_t u(const Vector &xvec)
{
   const int dim = xvec.Size();
   const real_t x = pi*xvec[0], y = pi*xvec[1];
   if (dim == 2) { return sin(x)*sin(y); }
   else { const real_t z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
}
 
std::function<real_t(const Vector &)> f(real_t mass_coeff)
{
   return [mass_coeff](const Vector &xvec)
   {
      const int dim = xvec.Size();
      const real_t x = pi*xvec[0], y = pi*xvec[1];
      if (dim == 2)
      {
         return mass_coeff*sin(x)*sin(y) + 2*pi2*sin(x)*sin(y);
      }
      else // dim == 3
      {
         const real_t z = pi*xvec[2];
         return mass_coeff*sin(x)*sin(y)*sin(z) + 3*pi2*sin(x)*sin(y)*sin(z);
      }
   };
}
 
// Exact solution for definite Maxwell and grad-div problems with RHS
// corresponding to f_vec below.
void u_vec(const Vector &xvec, Vector &u)
{
   const int dim = xvec.Size();
   const real_t x = pi*xvec[0], y = pi*xvec[1];
   if (dim == 2)
   {
      u[0] = cos(x)*sin(y);
      u[1] = sin(x)*cos(y);
   }
   else // dim == 3
   {
      const real_t z = pi*xvec[2];
      u[0] = cos(x)*sin(y)*sin(z);
      u[1] = sin(x)*cos(y)*sin(z);
      u[2] = sin(x)*sin(y)*cos(z);
   }
}
 
std::function<void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
{
   return [grad_div_problem](const Vector &xvec, Vector &f)
   {
      const int dim = xvec.Size();
      const real_t x = pi*xvec[0], y = pi*xvec[1];
      if (grad_div_problem)
      {
         if (dim == 2)
         {
            f[0] = (1 + 2*pi2)*cos(x)*sin(y);
            f[1] = (1 + 2*pi2)*cos(y)*sin(x);
         }
         else // dim == 3
         {
            const real_t z = pi*xvec[2];
            f[0] = (1 + 3*pi2)*cos(x)*sin(y)*sin(z);
            f[1] = (1 + 3*pi2)*cos(y)*sin(x)*sin(z);
            f[2] = (1 + 3*pi2)*cos(z)*sin(x)*sin(y);
         }
      }
      else
      {
         if (dim == 2)
         {
            f[0] = cos(x)*sin(y);
            f[1] = sin(x)*cos(y);
         }
         else // dim == 3
         {
            const real_t z = pi*xvec[2];
            f[0] = cos(x)*sin(y)*sin(z);
            f[1] = sin(x)*cos(y)*sin(z);
            f[2] = sin(x)*sin(y)*cos(z);
         }
      }
   };
}
 
} // namespace mfem
 
#endif