MFEM  v4.5.2
Finite element discretization library
mixed_integrator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_LIBCEED_MIXED_INTEGRATOR
13 #define MFEM_LIBCEED_MIXED_INTEGRATOR
14 
15 #include "ceed.hpp"
16 #include "integrator.hpp"
17 #include <unordered_map>
18 
19 namespace mfem
20 {
21 
22 namespace ceed
23 {
24 
25 /** @brief This class wraps a `ceed::PAIntegrator` or `ceed::MFIntegrator` to
26  support mixed finite element spaces. */
27 template <typename CeedInteg>
29 {
30 #ifdef MFEM_USE_CEED
31  using ElementKey = std::pair<int, int>; //< Element::Type, Order >
32  struct key_hash
33  {
34  std::size_t operator()(const ElementKey& k) const
35  {
36  return k.first + 2 * k.second;
37  }
38  };
39  using ElementsMap = std::unordered_map<const ElementKey, int*, key_hash>;
40  std::vector<CeedInteg*> sub_ops;
41 
42 public:
43  template <typename Integrator, typename CeedOperatorInfo, typename CoeffType>
44  void Assemble(const Integrator &integ,
45  CeedOperatorInfo &info,
46  const mfem::FiniteElementSpace &fes,
47  CoeffType *Q)
48  {
49  ElementsMap count;
50  ElementsMap element_indices;
51  ElementsMap offsets;
52 
53  // Count the number of elements of each type
54  for (int i = 0; i < fes.GetNE(); i++)
55  {
56  ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
57  auto value = count.find(key);
58  if (value == count.end())
59  {
60  count[key] = new int(1);
61  }
62  else
63  {
64  (*value->second)++;
65  }
66  }
67 
68  // Initialization of the arrays
69  for ( const auto& value : count )
70  {
71  element_indices[value.first] = new int[*value.second];
72  offsets[value.first] = new int(0);
73  }
74 
75  // Populates the indices arrays for each element type
76  for (int i = 0; i < fes.GetNE(); i++)
77  {
78  ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
79  int &offset = *(offsets[key]);
80  int* indices_array = element_indices[key];
81  indices_array[offset] = i;
82  offset++;
83  }
84 
85  // Create composite CeedOperator
86  CeedCompositeOperatorCreate(internal::ceed, &oper);
87 
88  // Create each sub-CeedOperator
89  sub_ops.reserve(element_indices.size());
90  for (const auto& value : element_indices)
91  {
92  const int* indices = value.second;
93  const int first_index = indices[0];
94  const mfem::FiniteElement &el = *fes.GetFE(first_index);
95  auto &T = *fes.GetMesh()->GetElementTransformation(first_index);
96  MFEM_ASSERT(!integ.GetIntegrationRule(),
97  "Mixed mesh integrators should not have an"
98  " IntegrationRule.");
99  const IntegrationRule &ir = GetRule(integ, el, el, T);
100  auto sub_op = new CeedInteg();
101  int nelem = *count[value.first];
102  sub_op->Assemble(info, fes, ir, nelem, indices, Q);
103  sub_ops.push_back(sub_op);
104  CeedCompositeOperatorAddSub(oper, sub_op->GetCeedOperator());
105  }
106 
107  const int ndofs = fes.GetVDim() * fes.GetNDofs();
108  CeedVectorCreate(internal::ceed, ndofs, &u);
109  CeedVectorCreate(internal::ceed, ndofs, &v);
110  }
111 
113  {
114  for (auto sub_op : sub_ops)
115  {
116  delete sub_op;
117  }
118  }
119 #endif
120 };
121 
122 } // namespace ceed
123 
124 } // namespace mfem
125 
126 #endif // MFEM_LIBCEED_MIXED_INTEGRATOR
Abstract class for all finite elements.
Definition: fe_base.hpp:232
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
int GetElementType(int i) const
Returns the type of element i.
Definition: fespace.hpp:635
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
CeedOperator oper
Definition: operator.hpp:29
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
This class wraps a ceed::PAIntegrator or ceed::MFIntegrator to support mixed finite element spaces...
void Assemble(const Integrator &integ, CeedOperatorInfo &info, const mfem::FiniteElementSpace &fes, CoeffType *Q)
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
const IntegrationRule & GetRule(const Integrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
int GetElementOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:178
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348