Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "basis.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "utils/diagnostic.hpp"
8 :
9 : namespace palace::ceed
10 : {
11 :
12 : namespace
13 : {
14 :
15 36650 : void InitTensorBasis(const mfem::FiniteElement &fe, const mfem::IntegrationRule &ir,
16 : CeedInt num_comp, Ceed ceed, CeedBasis *basis)
17 : {
18 : // The x-coordinates of the first `Q` points of the integration rule are the points of
19 : // the corresponding 1D rule. We also scale the weights accordingly.
20 36650 : const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
21 : const int dim = fe.GetDim();
22 36650 : const int P = maps.ndof;
23 36650 : const int Q = maps.nqpt;
24 : mfem::Vector qX(Q), qW(Q);
25 : double w_sum = 0.0;
26 158534 : for (int i = 0; i < Q; i++)
27 : {
28 : const mfem::IntegrationPoint &ip = ir.IntPoint(i);
29 121884 : qX(i) = ip.x;
30 121884 : qW(i) = ip.weight;
31 121884 : w_sum += ip.weight;
32 : }
33 36650 : qW *= 1.0 / w_sum;
34 :
35 36650 : PalaceCeedCall(ceed, CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, maps.Bt.GetData(),
36 : maps.Gt.GetData(), qX.GetData(),
37 : qW.GetData(), basis));
38 36650 : }
39 :
40 39551 : void InitNonTensorBasis(const mfem::FiniteElement &fe, const mfem::IntegrationRule &ir,
41 : CeedInt num_comp, Ceed ceed, CeedBasis *basis)
42 : {
43 39551 : const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::FULL);
44 : const int dim = fe.GetDim();
45 39551 : const int P = maps.ndof;
46 39551 : const int Q = maps.nqpt;
47 39551 : mfem::DenseMatrix qX(dim, Q);
48 : mfem::Vector qW(Q);
49 970182 : for (int i = 0; i < Q; i++)
50 : {
51 : const mfem::IntegrationPoint &ip = ir.IntPoint(i);
52 930631 : qX(0, i) = ip.x;
53 930631 : if (dim > 1)
54 : {
55 929659 : qX(1, i) = ip.y;
56 : }
57 929659 : if (dim > 2)
58 : {
59 728276 : qX(2, i) = ip.z;
60 : }
61 930631 : qW(i) = ip.weight;
62 : }
63 :
64 39551 : if (fe.GetMapType() == mfem::FiniteElement::H_DIV)
65 : {
66 4625 : PalaceCeedCall(ceed,
67 : CeedBasisCreateHdiv(ceed, GetCeedTopology(fe.GetGeomType()), num_comp, P,
68 : Q, maps.Bt.GetData(), maps.Gt.GetData(),
69 : qX.GetData(), qW.GetData(), basis));
70 : }
71 34926 : else if (fe.GetMapType() == mfem::FiniteElement::H_CURL)
72 : {
73 8273 : PalaceCeedCall(ceed,
74 : CeedBasisCreateHcurl(ceed, GetCeedTopology(fe.GetGeomType()), num_comp,
75 : P, Q, maps.Bt.GetData(), maps.Gt.GetData(),
76 : qX.GetData(), qW.GetData(), basis));
77 : }
78 : else
79 : {
80 26653 : PalaceCeedCall(ceed,
81 : CeedBasisCreateH1(ceed, GetCeedTopology(fe.GetGeomType()), num_comp, P,
82 : Q, maps.Bt.GetData(), maps.Gt.GetData(), qX.GetData(),
83 : qW.GetData(), basis));
84 : }
85 39551 : }
86 :
87 : PalacePragmaDiagnosticPush
88 : PalacePragmaDiagnosticDisableUnused
89 :
90 : void InitCeedInterpolatorBasis(const mfem::FiniteElement &trial_fe,
91 : const mfem::FiniteElement &test_fe, CeedInt trial_num_comp,
92 : CeedInt test_num_comp, Ceed ceed, CeedBasis *basis)
93 : {
94 : // Basis projection operator using libCEED.
95 : CeedBasis trial_basis, test_basis;
96 : const int P = std::max(trial_fe.GetDof(), test_fe.GetDof()), ir_order_max = 100;
97 : int ir_order = std::max(trial_fe.GetOrder(), test_fe.GetOrder());
98 : for (; ir_order < ir_order_max; ir_order++)
99 : {
100 : if (mfem::IntRules.Get(trial_fe.GetGeomType(), ir_order).GetNPoints() >= P)
101 : {
102 : break;
103 : }
104 : }
105 : const mfem::IntegrationRule &ir = mfem::IntRules.Get(trial_fe.GetGeomType(), ir_order);
106 :
107 : InitBasis(trial_fe, ir, trial_num_comp, ceed, &trial_basis);
108 : InitBasis(test_fe, ir, test_num_comp, ceed, &test_basis);
109 : PalaceCeedCall(ceed, CeedBasisCreateProjection(trial_basis, test_basis, basis));
110 : PalaceCeedCall(ceed, CeedBasisDestroy(&trial_basis));
111 : PalaceCeedCall(ceed, CeedBasisDestroy(&test_basis));
112 : }
113 :
114 : PalacePragmaDiagnosticPop
115 :
116 1254 : void InitMfemInterpolatorBasis(const mfem::FiniteElement &trial_fe,
117 : const mfem::FiniteElement &test_fe, CeedInt trial_num_comp,
118 : CeedInt test_num_comp, Ceed ceed, CeedBasis *basis)
119 : {
120 1254 : MFEM_VERIFY(trial_num_comp == test_num_comp && trial_num_comp == 1,
121 : "libCEED discrete linear operator requires same vdim = 1 for trial and test "
122 : "FE spaces!");
123 : const int trial_P = trial_fe.GetDof();
124 : const int test_P = test_fe.GetDof();
125 1254 : mfem::DenseMatrix Bt, Gt(trial_P, test_P);
126 : mfem::Vector qX(test_P), qW(test_P);
127 1254 : mfem::IsoparametricTransformation dummy;
128 1254 : dummy.SetIdentityTransformation(trial_fe.GetGeomType());
129 1254 : if (trial_fe.GetMapType() == test_fe.GetMapType())
130 : {
131 : // Prolongation.
132 816 : test_fe.GetTransferMatrix(trial_fe, dummy, Bt);
133 : }
134 438 : else if (trial_fe.GetMapType() == mfem::FiniteElement::VALUE &&
135 : test_fe.GetMapType() == mfem::FiniteElement::H_CURL)
136 : {
137 : // Discrete gradient interpolator.
138 303 : test_fe.ProjectGrad(trial_fe, dummy, Bt);
139 : }
140 135 : else if (trial_fe.GetMapType() == mfem::FiniteElement::H_CURL &&
141 : test_fe.GetMapType() == mfem::FiniteElement::H_DIV)
142 : {
143 : // Discrete curl interpolator.
144 135 : test_fe.ProjectCurl(trial_fe, dummy, Bt);
145 : }
146 0 : else if (trial_fe.GetMapType() == mfem::FiniteElement::H_DIV &&
147 : test_fe.GetMapType() == mfem::FiniteElement::INTEGRAL)
148 : {
149 : // Discrete divergence interpolator.
150 0 : test_fe.ProjectDiv(trial_fe, dummy, Bt);
151 : }
152 : else
153 : {
154 0 : MFEM_ABORT("Unsupported trial/test FE spaces for libCEED discrete linear operator!");
155 : }
156 1254 : Bt.Transpose();
157 1254 : Gt = 0.0;
158 1254 : qX = 0.0;
159 1254 : qW = 0.0;
160 :
161 : // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
162 1254 : PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, CEED_TOPOLOGY_LINE, trial_num_comp, trial_P,
163 : test_P, Bt.GetData(), Gt.GetData(), qX.GetData(),
164 : qW.GetData(), basis));
165 2508 : }
166 :
167 : } // namespace
168 :
169 76201 : void InitBasis(const mfem::FiniteElement &fe, const mfem::IntegrationRule &ir,
170 : CeedInt num_comp, Ceed ceed, CeedBasis *basis)
171 : {
172 : if constexpr (false)
173 : {
174 : std::cout << "New basis (" << ceed << ", " << &fe << ", " << &ir << ")\n";
175 : }
176 76201 : const bool tensor = dynamic_cast<const mfem::TensorBasisElement *>(&fe) != nullptr;
177 : const bool vector = fe.GetRangeType() == mfem::FiniteElement::VECTOR;
178 76201 : if (tensor && !vector)
179 : {
180 36650 : InitTensorBasis(fe, ir, num_comp, ceed, basis);
181 : }
182 : else
183 : {
184 39551 : InitNonTensorBasis(fe, ir, num_comp, ceed, basis);
185 : }
186 76201 : }
187 :
188 1254 : void InitInterpolatorBasis(const mfem::FiniteElement &trial_fe,
189 : const mfem::FiniteElement &test_fe, CeedInt trial_num_comp,
190 : CeedInt test_num_comp, Ceed ceed, CeedBasis *basis)
191 : {
192 : if constexpr (false)
193 : {
194 : std::cout << "New interpolator basis (" << ceed << ", " << &trial_fe << ", " << &test_fe
195 : << ", " << (trial_fe.GetMapType() == test_fe.GetMapType()) << ")\n";
196 : }
197 : if constexpr (false)
198 : {
199 : if (trial_fe.GetMapType() == test_fe.GetMapType())
200 : {
201 : InitCeedInterpolatorBasis(trial_fe, test_fe, trial_num_comp, test_num_comp, ceed,
202 : basis);
203 : }
204 : else
205 : {
206 : InitMfemInterpolatorBasis(trial_fe, test_fe, trial_num_comp, test_num_comp, ceed,
207 : basis);
208 : }
209 : }
210 : else
211 : {
212 1254 : InitMfemInterpolatorBasis(trial_fe, test_fe, trial_num_comp, test_num_comp, ceed,
213 : basis);
214 : }
215 1254 : }
216 :
217 : } // namespace palace::ceed
|
