finite_element.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "tuple.hpp"
#include "tensor.hpp"
#include "geometry.hpp"
#include "polynomials.hpp"
 
namespace serac {
 
template <int q>
struct TensorProductQuadratureRule {
  tensor<double, q> weights1D;  
  tensor<double, q> points1D;   
 
  SERAC_HOST_DEVICE double weight(int ix, int iy) const { return weights1D[ix] * weights1D[iy]; }
 
  SERAC_HOST_DEVICE double weight(int ix, int iy, int iz) const
  {
    return weights1D[ix] * weights1D[iy] * weights1D[iz];
  }
};
 
template <auto val>
struct CompileTimeValue {
};
 
template <mfem::Geometry::Type g, int q>
struct batched_jacobian;
 
template <int q>
struct batched_jacobian<mfem::Geometry::CUBE, q> {
  using type = tensor<double, 3, 3, q * q * q>;
};
 
template <int q>
struct batched_jacobian<mfem::Geometry::SQUARE, q> {
  using type = tensor<double, 2, 2, q * q>;
};
 
template <int q>
struct batched_jacobian<mfem::Geometry::TRIANGLE, q> {
  using type = tensor<double, 2, 2, q*(q + 1) / 2>;
};
 
template <int q>
struct batched_jacobian<mfem::Geometry::TETRAHEDRON, q> {
  using type = tensor<double, 3, 3, (q * (q + 1) * (q + 2)) / 6>;
};
 
template <mfem::Geometry::Type g, int q>
struct batched_position;
 
template <int q>
struct batched_position<mfem::Geometry::CUBE, q> {
  using type = tensor<double, 3, q * q * q>;
};
 
template <int q>
struct batched_position<mfem::Geometry::SQUARE, q> {
  using type = tensor<double, 2, q * q>;
};
 
template <int q>
struct batched_position<mfem::Geometry::TRIANGLE, q> {
  using type = tensor<double, 2, q*(q + 1) / 2>;
};
 
template <int q>
struct batched_position<mfem::Geometry::TETRAHEDRON, q> {
  using type = tensor<double, 3, (q * (q + 1) * (q + 2)) / 6>;
};
 
template <int q>
struct batched_position<mfem::Geometry::SEGMENT, q> {
  using type = tensor<double, q>;
};
 
template <mfem::Geometry::Type g>
SERAC_HOST_DEVICE constexpr int elements_per_block(int q)
{
  if (g == mfem::Geometry::CUBE) {
    switch (q) {
      case 1:
        return 64;
      case 2:
        return 16;
      case 3:
        return 4;
      default:
        return 1;
    }
  }
 
  if (g == mfem::Geometry::SQUARE) {
    switch (q) {
      case 1:
        return 128;
      case 2:
        return 32;
      case 3:
        return 16;
      case 4:
        return 8;
      default:
        return 1;
    }
  }
}
 
enum class Family
{
  QOI,
  H1,
  HCURL,
  HDIV,
  L2
};
 
template <int p, int c = 1>
struct H1 {
  static constexpr int    order      = p;           
  static constexpr int    components = c;           
  static constexpr Family family     = Family::H1;  
};
 
template <int p, int c = 1>
struct Hcurl {
  static constexpr int    order      = p;              
  static constexpr int    components = c;              
  static constexpr Family family     = Family::HCURL;  
};
 
template <int p, int c = 1>
struct L2 {
  static constexpr int    order      = p;           
  static constexpr int    components = c;           
  static constexpr Family family     = Family::L2;  
};
 
struct QOI {
  static constexpr int    order      = 0;            
  static constexpr int    components = 1;            
  static constexpr Family family     = Family::QOI;  
};
 
template <Family f, typename T, int q, int dim>
SERAC_HOST_DEVICE void parent_to_physical(tensor<T, q>& qf_input, const tensor<double, dim, dim, q>& jacobians)
{
  [[maybe_unused]] constexpr int VALUE      = 0;
  [[maybe_unused]] constexpr int DERIVATIVE = 1;
 
  for (int k = 0; k < q; k++) {
    tensor<double, dim, dim> J;
    for (int row = 0; row < dim; row++) {
      for (int col = 0; col < dim; col++) {
        J[row][col] = jacobians(col, row, k);
      }
    }
 
    if constexpr (f == Family::H1 || f == Family::L2) {
      // note: no transformation necessary for the values of H1-field
      get<DERIVATIVE>(qf_input[k]) = dot(get<DERIVATIVE>(qf_input[k]), inv(J));
    }
 
    if constexpr (f == Family::HCURL) {
      get<VALUE>(qf_input[k])      = dot(get<VALUE>(qf_input[k]), inv(J));
      get<DERIVATIVE>(qf_input[k]) = get<DERIVATIVE>(qf_input[k]) / det(J);
      if constexpr (dim == 3) {
        get<DERIVATIVE>(qf_input[k]) = dot(get<DERIVATIVE>(qf_input[k]), transpose(J));
      }
    }
  }
}
 
template <Family f, typename T, int q, int dim>
SERAC_HOST_DEVICE void physical_to_parent(tensor<T, q>& qf_output, const tensor<double, dim, dim, q>& jacobians)
{
  [[maybe_unused]] constexpr int SOURCE = 0;
  [[maybe_unused]] constexpr int FLUX   = 1;
 
  for (int k = 0; k < q; k++) {
    tensor<double, dim, dim> J_T;
    for (int row = 0; row < dim; row++) {
      for (int col = 0; col < dim; col++) {
        J_T[row][col] = jacobians(row, col, k);
      }
    }
 
    auto dv = det(J_T);
 
    if constexpr (f == Family::H1 || f == Family::L2) {
      get<SOURCE>(qf_output[k]) = get<SOURCE>(qf_output[k]) * dv;
      get<FLUX>(qf_output[k])   = dot(get<FLUX>(qf_output[k]), inv(J_T)) * dv;
    }
 
    // note: the flux term here is usually divided by detJ, but
    // physical_to_parent also multiplies every quadrature-point value by det(J)
    // so that part cancels out
    if constexpr (f == Family::HCURL) {
      get<SOURCE>(qf_output[k]) = dot(get<SOURCE>(qf_output[k]), inv(J_T)) * dv;
      if constexpr (dim == 3) {
        get<FLUX>(qf_output[k]) = dot(get<FLUX>(qf_output[k]), transpose(J_T));
      }
    }
 
    if constexpr (f == Family::QOI) {
      qf_output[k] = qf_output[k] * dv;
    }
  }
}
 
template <mfem::Geometry::Type g, typename family>
struct finite_element;
 
#include "detail/segment_H1.inl"
#include "detail/segment_Hcurl.inl"
#include "detail/segment_L2.inl"
 
#include "detail/triangle_H1.inl"
#include "detail/triangle_L2.inl"
 
#include "detail/quadrilateral_H1.inl"
#include "detail/quadrilateral_Hcurl.inl"
#include "detail/quadrilateral_L2.inl"
 
#include "detail/tetrahedron_H1.inl"
#include "detail/tetrahedron_L2.inl"
 
#include "detail/hexahedron_H1.inl"
#include "detail/hexahedron_Hcurl.inl"
#include "detail/hexahedron_L2.inl"
 
#include "detail/qoi.inl"
 
}  // namespace serac