domain.hpp

Go to the documentation of this file.

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <cstdlib>
#include <vector>
#include <array>
#include <cstdint>
#include <functional>
#include <map>
#include <set>
 
#include "mfem.hpp"
 
#include "smith/numerics/functional/tensor.hpp"
#include "smith/numerics/functional/finite_element.hpp"
#include "smith/numerics/functional/element_restriction.hpp"
#include "smith/numerics/functional/typedefs.hpp"
 
namespace smith {
 
struct BlockElementRestriction;
 
struct Domain {
  enum Type
  {
    Elements,
    BoundaryElements,
    InteriorBoundaryElements
  };
 
  static constexpr int num_types = 3;  
 
  const mesh_t& mesh_;
 
  int dim_;
 
  Type type_;
 
 
  std::vector<int> edge_ids_;
  std::vector<int> tri_ids_;
  std::vector<int> quad_ids_;
  std::vector<int> tet_ids_;
  std::vector<int> hex_ids_;
 
  std::vector<int> mfem_edge_ids_;
  std::vector<int> mfem_tri_ids_;
  std::vector<int> mfem_quad_ids_;
  std::vector<int> mfem_tet_ids_;
  std::vector<int> mfem_hex_ids_;
 
  std::map<FunctionSpace, BlockElementRestriction> restriction_operators;
 
  std::vector<int> shared_interior_face_ids_;
 
  Domain(const mesh_t& m, int d, Type type) : mesh_(m), dim_(d), type_(type) {}
 
  static Domain ofVertices(const mesh_t& mesh, std::function<bool(vec2)> func);
 
  static Domain ofVertices(const mesh_t& mesh, std::function<bool(vec3)> func);
 
  static Domain ofEdges(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
 
  static Domain ofEdges(const mesh_t& mesh, std::function<bool(std::vector<vec3>)> func);
 
  static Domain ofFaces(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
 
  static Domain ofFaces(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
 
  static Domain ofElements(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
 
  static Domain ofElements(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
 
  static Domain ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
 
  static Domain ofBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
 
  static Domain ofInteriorBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec2>, int)> func);
 
  static Domain ofInteriorBoundaryElements(const mesh_t& mesh, std::function<bool(std::vector<vec3>, int)> func);
 
  const std::vector<int>& get(mfem::Geometry::Type geom) const
  {
    if (geom == mfem::Geometry::SEGMENT) return edge_ids_;
    if (geom == mfem::Geometry::TRIANGLE) return tri_ids_;
    if (geom == mfem::Geometry::SQUARE) return quad_ids_;
    if (geom == mfem::Geometry::TETRAHEDRON) return tet_ids_;
    if (geom == mfem::Geometry::CUBE) return hex_ids_;
 
    exit(1);
  }
 
  const std::vector<int>& get_mfem_ids(mfem::Geometry::Type geom) const
  {
    if (geom == mfem::Geometry::SEGMENT) return mfem_edge_ids_;
    if (geom == mfem::Geometry::TRIANGLE) return mfem_tri_ids_;
    if (geom == mfem::Geometry::SQUARE) return mfem_quad_ids_;
    if (geom == mfem::Geometry::TETRAHEDRON) return mfem_tet_ids_;
    if (geom == mfem::Geometry::CUBE) return mfem_hex_ids_;
 
    exit(1);
  }
 
  int total_elements() const
  {
    return int(edge_ids_.size() + tri_ids_.size() + quad_ids_.size() + tet_ids_.size() + hex_ids_.size());
  }
 
  mfem::Array<int> bOffsets() const
  {
    mfem::Array<int> offsets(mfem::Geometry::NUM_GEOMETRIES + 1);
 
    int total = 0;
    offsets[mfem::Geometry::POINT] = total;
    total += 0;  // vertices;
    offsets[mfem::Geometry::SEGMENT] = total;
    total += int(edge_ids_.size());
    offsets[mfem::Geometry::TRIANGLE] = total;
    total += int(tri_ids_.size());
    offsets[mfem::Geometry::SQUARE] = total;
    total += int(quad_ids_.size());
    offsets[mfem::Geometry::TETRAHEDRON] = total;
    total += int(tet_ids_.size());
    offsets[mfem::Geometry::CUBE] = total;
    total += int(hex_ids_.size());
    offsets[mfem::Geometry::PRISM] = total;
    offsets[mfem::Geometry::PYRAMID] = total;
    offsets[mfem::Geometry::NUM_GEOMETRIES] = total;
 
    return offsets;
  }
 
  mfem::Array<int> dof_list(const fes_t* fes) const;
 
  void insert_restriction(const fes_t* fes, FunctionSpace space);
 
  const BlockElementRestriction& get_restriction(FunctionSpace space);
 
  void addElement(int geom_id, int elem_id, mfem::Geometry::Type element_geometry);
 
  void addElements(const std::vector<int>& geom_id, const std::vector<int>& elem_id,
                   mfem::Geometry::Type element_geometry);
 
  void insert_shared_interior_face_list();
};
 
Domain EntireDomain(const mesh_t& mesh);
 
Domain EntireBoundary(const mesh_t& mesh);
 
Domain EntireInteriorBoundary(const mesh_t& mesh);
 
Domain operator|(const Domain& a, const Domain& b);
 
Domain operator&(const Domain& a, const Domain& b);
 
Domain operator-(const Domain& a, const Domain& b);
 
template <int dim>
inline auto by_attr(int value)
{
  return [value](std::vector<tensor<double, dim>>, int attr) { return value == attr; };
}
 
template <int dim>
inline auto by_attr(std::set<int> values)
{
  return [values](std::vector<tensor<double, dim>>, int attr) { return values.find(attr) != values.end(); };
}
 
inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(const Domain& domain)
{
  std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> counts{};
 
  constexpr std::array<mfem::Geometry::Type, 5> geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE,
                                                              mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON,
                                                              mfem::Geometry::CUBE};
  for (auto geom : geometries) {
    counts[uint32_t(geom)] = uint32_t(domain.get(geom).size());
  }
  return counts;
}
 
template <int dim>
inline tensor<double, dim> average(std::vector<tensor<double, dim>>& positions)
{
  tensor<double, dim> total{};
  for (auto x : positions) {
    total += x;
  }
  return total / double(positions.size());
}
 
}  // namespace smith