MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
basis.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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 "../../gridfunc.hpp"
13#include "util.hpp"
14
15namespace mfem
16{
17
18namespace ceed
19{
20
21#ifdef MFEM_USE_CEED
22
23static CeedElemTopology GetCeedTopology(Geometry::Type geom)
24{
25 switch (geom)
26 {
28 return CEED_TOPOLOGY_LINE;
30 return CEED_TOPOLOGY_TRIANGLE;
32 return CEED_TOPOLOGY_QUAD;
34 return CEED_TOPOLOGY_TET;
35 case Geometry::CUBE:
36 return CEED_TOPOLOGY_HEX;
37 case Geometry::PRISM:
38 return CEED_TOPOLOGY_PRISM;
40 return CEED_TOPOLOGY_PYRAMID;
41 default:
42 MFEM_ABORT("This type of element is not supported");
43 return CEED_TOPOLOGY_PRISM; // Silence warning
44 }
45}
46
47static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
48 const mfem::FiniteElement &fe,
49 const mfem::IntegrationRule &ir,
50 Ceed ceed, CeedBasis *basis)
51{
53 mfem::Mesh *mesh = fes.GetMesh();
54 const int dim = mesh->Dimension();
55 const int ndofs = maps.ndof;
56 const int nqpts = maps.nqpt;
57 mfem::DenseMatrix qX(dim,nqpts);
58 mfem::Vector qW(nqpts);
59 for (int i = 0; i < nqpts; i++)
60 {
61 const mfem::IntegrationPoint &ip = ir.IntPoint(i);
62 qX(0,i) = ip.x;
63 if (dim>1) { qX(1,i) = ip.y; }
64 if (dim>2) { qX(2,i) = ip.z; }
65 qW(i) = ip.weight;
66 }
67 CeedBasisCreateH1(ceed, GetCeedTopology(fe.GetGeomType()),
68 fes.GetVDim(), ndofs, nqpts,
69 maps.Bt.GetData(), maps.Gt.GetData(),
70 qX.GetData(), qW.GetData(), basis);
71}
72
73static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
74 const mfem::FiniteElement &fe,
75 const mfem::IntegrationRule &ir,
76 Ceed ceed, CeedBasis *basis)
77{
79 mfem::Mesh *mesh = fes.GetMesh();
80 const int ndofs = maps.ndof;
81 const int nqpts = maps.nqpt;
82 mfem::Vector qX(nqpts), qW(nqpts);
83 // The x-coordinates of the first `nqpts` points of the integration rule are
84 // the points of the corresponding 1D rule. We also scale the weights
85 // accordingly.
86 double w_sum = 0.0;
87 for (int i = 0; i < nqpts; i++)
88 {
89 const mfem::IntegrationPoint &ip = ir.IntPoint(i);
90 qX(i) = ip.x;
91 qW(i) = ip.weight;
92 w_sum += ip.weight;
93 }
94 qW *= 1.0/w_sum;
95 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes.GetVDim(), ndofs,
96 nqpts, maps.Bt.GetData(),
97 maps.Gt.GetData(), qX.GetData(),
98 qW.GetData(), basis);
99}
100
101static void InitBasisImpl(const FiniteElementSpace &fes,
102 const FiniteElement &fe,
103 const IntegrationRule &ir,
104 Ceed ceed, CeedBasis *basis)
105{
106 // Check for FES -> basis, restriction in hash tables
107 const int P = fe.GetDof();
108 const int Q = ir.GetNPoints();
109 const int ncomp = fes.GetVDim();
110 BasisKey basis_key(&fes, &ir, ncomp, P, Q);
111 auto basis_itr = mfem::internal::ceed_basis_map.find(basis_key);
112 const bool tensor = dynamic_cast<const mfem::TensorBasisElement *>
113 (&fe) != nullptr;
114
115 // Init or retrieve key values
116 if (basis_itr == mfem::internal::ceed_basis_map.end())
117 {
118 if ( tensor )
119 {
120 InitTensorBasis(fes, fe, ir, ceed, basis);
121 }
122 else
123 {
124 InitNonTensorBasis(fes, fe, ir, ceed, basis);
125 }
126 mfem::internal::ceed_basis_map[basis_key] = *basis;
127 }
128 else
129 {
130 *basis = basis_itr->second;
131 }
132}
133
135 const IntegrationRule &ir,
136 Ceed ceed, CeedBasis *basis)
137{
138 const mfem::FiniteElement &fe = *fes.GetFE(0);
139 InitBasisImpl(fes, fe, ir, ceed, basis);
140}
141
143 const IntegrationRule &ir,
144 int nelem,
145 const int* indices,
146 Ceed ceed, CeedBasis *basis)
147{
148 const mfem::FiniteElement &fe = *fes.GetFE(indices[0]);
149 InitBasisImpl(fes, fe, ir, ceed, basis);
150}
151
152#endif
153
154} // namespace ceed
155
156} // namespace mfem
T * GetData()
Returns the data.
Definition array.hpp:118
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:137
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:154
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:161
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:174
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:217
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:195
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Abstract class for all finite elements.
Definition fe_base.hpp:239
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:326
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Mesh data type.
Definition mesh.hpp:56
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
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
std::tuple< const mfem::FiniteElementSpace *, const mfem::IntegrationRule *, int, int, int > BasisKey
Definition util.hpp:127
void InitBasis(const FiniteElementSpace &fes, const IntegrationRule &ir, Ceed ceed, CeedBasis *basis)
Initialize a CeedBasis for non-mixed meshes.
Definition basis.cpp:134