MFEM  v4.6.0
Finite element discretization library
util.cpp
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 #include "../../../general/device.hpp"
13 #include "../../../fem/gridfunc.hpp"
14 #include "../../../linalg/dtensor.hpp"
15 
16 #include "basis.hpp"
17 #include "restriction.hpp"
18 #include "ceed.hpp"
19 
20 #include <sys/types.h>
21 #include <sys/stat.h>
22 #if !defined(_WIN32) || !defined(_MSC_VER)
23 typedef struct stat struct_stat;
24 #else
25 #define stat(dir, buf) _stat(dir, buf)
26 #define S_ISDIR(mode) (((mode) & _S_IFMT) == _S_IFDIR)
27 typedef struct _stat struct_stat;
28 #endif
29 
30 namespace mfem
31 {
32 
34 {
36 }
37 
38 namespace ceed
39 {
40 
42 {
43 #ifdef MFEM_USE_CEED
44  auto itb = mfem::internal::ceed_basis_map.begin();
45  while (itb != mfem::internal::ceed_basis_map.end())
46  {
47  if (std::get<0>(itb->first)==fes)
48  {
49  CeedBasisDestroy(&itb->second);
50  itb = mfem::internal::ceed_basis_map.erase(itb);
51  }
52  else
53  {
54  itb++;
55  }
56  }
57  auto itr = mfem::internal::ceed_restr_map.begin();
58  while (itr != mfem::internal::ceed_restr_map.end())
59  {
60  if (std::get<0>(itr->first)==fes)
61  {
62  CeedElemRestrictionDestroy(&itr->second);
63  itr = mfem::internal::ceed_restr_map.erase(itr);
64  }
65  else
66  {
67  itr++;
68  }
69  }
70 #endif
71 }
72 
73 #ifdef MFEM_USE_CEED
74 
75 void InitVector(const mfem::Vector &v, CeedVector &cv)
76 {
77  CeedVectorCreate(mfem::internal::ceed, v.Size(), &cv);
78  CeedScalar *cv_ptr;
79  CeedMemType mem;
80  CeedGetPreferredMemType(mfem::internal::ceed, &mem);
81  if ( Device::Allows(Backend::DEVICE_MASK) && mem==CEED_MEM_DEVICE )
82  {
83  cv_ptr = const_cast<CeedScalar*>(v.Read());
84  }
85  else
86  {
87  cv_ptr = const_cast<CeedScalar*>(v.HostRead());
88  mem = CEED_MEM_HOST;
89  }
90  CeedVectorSetArray(cv, mem, CEED_USE_POINTER, cv_ptr);
91 }
92 
94  const IntegrationRule &irm,
95  Ceed ceed, CeedBasis *basis,
96  CeedElemRestriction *restr)
97 {
98  InitBasis(fes, irm, ceed, basis);
99  InitRestriction(fes, ceed, restr);
100 }
101 
103  const IntegrationRule &irm,
104  int nelem,
105  const int* indices,
106  Ceed ceed, CeedBasis *basis,
107  CeedElemRestriction *restr)
108 {
109  InitBasisWithIndices(fes, irm, nelem, indices, ceed, basis);
110  InitRestrictionWithIndices(fes, nelem, indices, ceed, restr);
111 }
112 
114  const IntegrationRule &irm,
115  int nelem,
116  const int* indices,
117  Ceed ceed, CeedBasis *basis,
118  CeedElemRestriction *restr)
119 {
120  if (indices)
121  {
122  InitBasisAndRestrictionWithIndices(fes,irm,nelem,indices,ceed,basis,restr);
123  }
124  else
125  {
126  InitBasisAndRestriction(fes,irm,ceed,basis,restr);
127  }
128 }
129 
130 // Assumes a tensor-product operator with one active field
131 int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
132 {
133  int ierr;
134  Ceed ceed;
135  ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
136 
137  CeedQFunction qf;
138  bool isComposite;
139  ierr = CeedOperatorIsComposite(oper, &isComposite); CeedChk(ierr);
140  CeedOperator *subops;
141  if (isComposite)
142  {
143 #if CEED_VERSION_GE(0, 10, 2)
144  ierr = CeedCompositeOperatorGetSubList(oper, &subops); CeedChk(ierr);
145 #else
146  ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
147 #endif
148  ierr = CeedOperatorGetQFunction(subops[0], &qf); CeedChk(ierr);
149  }
150  else
151  {
152  ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
153  }
154  CeedInt numinputfields, numoutputfields;
155  ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
156  CeedOperatorField *inputfields;
157  if (isComposite)
158  {
159  ierr = CeedOperatorGetFields(subops[0], &numinputfields, &inputfields,
160  &numoutputfields, NULL); CeedChk(ierr);
161  }
162  else
163  {
164  ierr = CeedOperatorGetFields(oper, &numinputfields, &inputfields,
165  &numoutputfields, NULL); CeedChk(ierr);
166  }
167 
168  CeedVector if_vector;
169  bool found = false;
170  int found_index = -1;
171  for (int i = 0; i < numinputfields; ++i)
172  {
173  ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector); CeedChk(ierr);
174  if (if_vector == CEED_VECTOR_ACTIVE)
175  {
176  if (found)
177  {
178  return CeedError(ceed, 1, "Multiple active vectors in CeedOperator!");
179  }
180  found = true;
181  found_index = i;
182  }
183  }
184  if (!found)
185  {
186  return CeedError(ceed, 1, "No active vector in CeedOperator!");
187  }
188  *field = inputfields[found_index];
189 
190  return 0;
191 }
192 
193 template <>
195  const MassIntegrator &integ,
196  const FiniteElement &trial_fe,
197  const FiniteElement &test_fe,
199 {
200  return MassIntegrator::GetRule(trial_fe, test_fe, trans);
201 }
202 
203 template <>
205  const VectorMassIntegrator &integ,
206  const FiniteElement &trial_fe,
207  const FiniteElement &test_fe,
209 {
210  return MassIntegrator::GetRule(trial_fe, test_fe, trans);
211 }
212 
213 template <>
215  const ConvectionIntegrator &integ,
216  const FiniteElement &trial_fe,
217  const FiniteElement &test_fe,
219 {
220  return ConvectionIntegrator::GetRule(trial_fe, test_fe, trans);
221 }
222 
223 template <>
225  const VectorConvectionNLFIntegrator &integ,
226  const FiniteElement &trial_fe,
227  const FiniteElement &test_fe,
229 {
231 }
232 
233 template <>
235  const DiffusionIntegrator &integ,
236  const FiniteElement &trial_fe,
237  const FiniteElement &test_fe,
239 {
240  return DiffusionIntegrator::GetRule(trial_fe, test_fe);
241 }
242 
243 template <>
245  const VectorDiffusionIntegrator &integ,
246  const FiniteElement &trial_fe,
247  const FiniteElement &test_fe,
249 {
250  return DiffusionIntegrator::GetRule(trial_fe, test_fe);
251 }
252 
253 std::string ceed_path;
254 
255 const std::string &GetCeedPath()
256 {
257  if (ceed_path.empty())
258  {
259  const char *install_dir = MFEM_INSTALL_DIR "/include/mfem/fem/ceed";
260  const char *source_dir = MFEM_SOURCE_DIR "/fem/ceed";
261  struct_stat m_stat;
262  if (stat(install_dir, &m_stat) == 0 && S_ISDIR(m_stat.st_mode))
263  {
264  ceed_path = install_dir;
265  }
266  else if (stat(source_dir, &m_stat) == 0 && S_ISDIR(m_stat.st_mode))
267  {
268  ceed_path = source_dir;
269  }
270  else
271  {
272  MFEM_ABORT("Cannot find libCEED kernels in MFEM_INSTALL_DIR or "
273  "MFEM_SOURCE_DIR");
274  }
275  // Could be useful for debugging:
276  // out << "Using libCEED dir: " << ceed_path << std::endl;
277  }
278  return ceed_path;
279 }
280 
281 #endif
282 
283 } // namespace ceed
284 
285 } // namespace mfem
void InitBasisWithIndices(const FiniteElementSpace &fes, const IntegrationRule &ir, int nelem, const int *indices, Ceed ceed, CeedBasis *basis)
Initialize a CeedBasis for mixed meshes.
Definition: basis.cpp:142
Abstract class for all finite elements.
Definition: fe_base.hpp:233
void InitBasis(const FiniteElementSpace &fes, const IntegrationRule &ir, Ceed ceed, CeedBasis *basis)
Initialize a CeedBasis for non-mixed meshes.
Definition: basis.cpp:134
std::string ceed_path
Definition: util.cpp:253
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
const IntegrationRule & GetRule< ConvectionIntegrator >(const ConvectionIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:214
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
void InitVector(const mfem::Vector &v, CeedVector &cv)
Initialize a CeedVector from an mfem::Vector.
Definition: util.cpp:75
struct stat struct_stat
Definition: util.cpp:23
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
void InitBasisAndRestrictionWithIndices(const FiniteElementSpace &fes, const IntegrationRule &irm, int nelem, const int *indices, Ceed ceed, CeedBasis *basis, CeedElemRestriction *restr)
Definition: util.cpp:102
int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
Definition: util.cpp:131
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
const IntegrationRule & GetRule< VectorDiffusionIntegrator >(const VectorDiffusionIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:244
const std::string & GetCeedPath()
Return the path to the libCEED q-function headers.
Definition: util.cpp:255
static const IntegrationRule & GetRule(const FiniteElement &el, ElementTransformation &Trans)
void RemoveBasisAndRestriction(const mfem::FiniteElementSpace *fes)
Remove from ceed_basis_map and ceed_restr_map the entries associated with the given fes...
Definition: util.cpp:41
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
Definition: util.cpp:33
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe)
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Definition: device.hpp:258
void InitRestrictionWithIndices(const FiniteElementSpace &fes, int nelem, const int *indices, Ceed ceed, CeedElemRestriction *restr)
Initialize a CeedElemRestriction for mixed meshes.
Vector data type.
Definition: vector.hpp:58
const IntegrationRule & GetRule< VectorConvectionNLFIntegrator >(const VectorConvectionNLFIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:224
const IntegrationRule & GetRule< DiffusionIntegrator >(const DiffusionIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:234
void InitBasisAndRestriction(const FiniteElementSpace &fes, const IntegrationRule &irm, Ceed ceed, CeedBasis *basis, CeedElemRestriction *restr)
Initialize a CeedBasis and a CeedElemRestriction based on an mfem::FiniteElementSpace fes...
Definition: util.cpp:93
Biwise-OR of all device backends.
Definition: device.hpp:96
static const IntegrationRule & GetRule(const FiniteElement &fe, ElementTransformation &T)
const IntegrationRule & GetRule< MassIntegrator >(const MassIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:194
Bitwise-OR of all CEED backends.
Definition: device.hpp:94
void InitRestriction(const FiniteElementSpace &fes, Ceed ceed, CeedElemRestriction *restr)
Initialize a CeedElemRestriction for non-mixed meshes.
const IntegrationRule & GetRule< VectorMassIntegrator >(const VectorMassIntegrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &trans)
Definition: util.cpp:204
alpha (q . grad u, v)