MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
full-assembly.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 "full-assembly.hpp"
13
14#ifdef MFEM_USE_CEED
15
16#include "../../../linalg/sparsemat.hpp"
17#include "../interface/util.hpp"
18#include "../interface/ceed.hpp"
19
20namespace mfem
21{
22
23namespace ceed
24{
25
26int CeedHackReallocArray(size_t n, size_t unit, void *p)
27{
28 *(void **)p = realloc(*(void **)p, n*unit);
29 if (n && unit && !*(void **)p)
30 return CeedError(NULL, 1, "realloc failed to allocate %zd members of size "
31 "%zd\n", n, unit);
32 return 0;
33}
34
35#define CeedHackRealloc(n, p) CeedHackReallocArray((n), sizeof(**(p)), p)
36
37int CeedHackFree(void *p)
38{
39 free(*(void **)p);
40 *(void **)p = NULL;
41 return 0;
42}
43
44int CeedSingleOperatorFullAssemble(CeedOperator op, SparseMatrix *out)
45{
46 int ierr;
47 Ceed ceed;
48 ierr = CeedOperatorGetCeed(op, &ceed); PCeedChk(ierr);
49
50 // Assemble QFunction
51 CeedQFunction qf;
52 ierr = CeedOperatorGetQFunction(op, &qf); PCeedChk(ierr);
53 CeedInt numinputfields, numoutputfields;
54 PCeedChk(ierr);
55 CeedVector assembledqf;
56 CeedElemRestriction rstr_q;
57 ierr = CeedOperatorLinearAssembleQFunction(
58 op, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); PCeedChk(ierr);
59
60 CeedSize qflength;
61 ierr = CeedVectorGetLength(assembledqf, &qflength); PCeedChk(ierr);
62
63 CeedOperatorField *input_fields;
64 CeedOperatorField *output_fields;
65 ierr = CeedOperatorGetFields(op, &numinputfields, &input_fields,
66 &numoutputfields, &output_fields);
67 PCeedChk(ierr);
68
69 // Determine active input basis
70 CeedQFunctionField *qffields;
71 ierr = CeedQFunctionGetFields(qf, &numinputfields, &qffields,
72 &numoutputfields, NULL);
73 PCeedChk(ierr);
74 CeedInt numemodein = 0, ncomp, dim = 1;
75 CeedEvalMode *emodein = NULL;
76 CeedBasis basisin = NULL;
77 CeedElemRestriction rstrin = NULL;
78 for (CeedInt i=0; i<numinputfields; i++)
79 {
80 CeedVector vec;
81 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); PCeedChk(ierr);
82 if (vec == CEED_VECTOR_ACTIVE)
83 {
84 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basisin);
85 PCeedChk(ierr);
86 ierr = CeedBasisGetNumComponents(basisin, &ncomp); PCeedChk(ierr);
87 ierr = CeedBasisGetDimension(basisin, &dim); PCeedChk(ierr);
88 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstrin);
89 PCeedChk(ierr);
90 CeedEvalMode emode;
91 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
92 PCeedChk(ierr);
93 switch (emode)
94 {
95 case CEED_EVAL_NONE:
96 case CEED_EVAL_INTERP:
97 ierr = CeedHackRealloc(numemodein + 1, &emodein); PCeedChk(ierr);
98 emodein[numemodein] = emode;
99 numemodein += 1;
100 break;
101 case CEED_EVAL_GRAD:
102 ierr = CeedHackRealloc(numemodein + dim, &emodein); PCeedChk(ierr);
103 for (CeedInt d=0; d<dim; d++)
104 {
105 emodein[numemodein+d] = emode;
106 }
107 numemodein += dim;
108 break;
109 case CEED_EVAL_WEIGHT:
110 case CEED_EVAL_DIV:
111 case CEED_EVAL_CURL:
112 break; // Caught by QF Assembly
113 }
114 }
115 }
116
117 // Determine active output basis
118 ierr = CeedQFunctionGetFields(qf, &numinputfields, NULL, &numoutputfields,
119 &qffields); PCeedChk(ierr);
120 CeedInt numemodeout = 0;
121 CeedEvalMode *emodeout = NULL;
122 CeedBasis basisout = NULL;
123 CeedElemRestriction rstrout = NULL;
124 for (CeedInt i=0; i<numoutputfields; i++)
125 {
126 CeedVector vec;
127 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); PCeedChk(ierr);
128 if (vec == CEED_VECTOR_ACTIVE)
129 {
130 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basisout);
131 PCeedChk(ierr);
132 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstrout);
133 PCeedChk(ierr);
134 PCeedChk(ierr);
135 CeedEvalMode emode;
136 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
137 PCeedChk(ierr);
138 switch (emode)
139 {
140 case CEED_EVAL_NONE:
141 case CEED_EVAL_INTERP:
142 ierr = CeedHackRealloc(numemodeout + 1, &emodeout); PCeedChk(ierr);
143 emodeout[numemodeout] = emode;
144 numemodeout += 1;
145 break;
146 case CEED_EVAL_GRAD:
147 ierr = CeedHackRealloc(numemodeout + dim, &emodeout); PCeedChk(ierr);
148 for (CeedInt d=0; d<dim; d++)
149 {
150 emodeout[numemodeout+d] = emode;
151 }
152 numemodeout += dim;
153 break;
154 case CEED_EVAL_WEIGHT:
155 case CEED_EVAL_DIV:
156 case CEED_EVAL_CURL:
157 break; // Caught by QF Assembly
158 }
159 }
160 }
161
162 CeedInt nelem, elemsize, nqpts;
163 CeedSize nnodes;
164 ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); PCeedChk(ierr);
165 ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); PCeedChk(ierr);
166 ierr = CeedElemRestrictionGetLVectorSize(rstrin, &nnodes); PCeedChk(ierr);
167 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); PCeedChk(ierr);
168
169 // Determine elem_dof relation
170 CeedVector index_vec;
171 ierr = CeedVectorCreate(ceed, nnodes, &index_vec); PCeedChk(ierr);
172 CeedScalar *array;
173 ierr = CeedVectorGetArrayWrite(index_vec, CEED_MEM_HOST, &array);
174 PCeedChk(ierr);
175 for (CeedSize i = 0; i < nnodes; ++i)
176 {
177 array[i] = i;
178 }
179 ierr = CeedVectorRestoreArray(index_vec, &array); PCeedChk(ierr);
180 CeedVector elem_dof;
181 ierr = CeedVectorCreate(ceed, nelem * elemsize, &elem_dof); PCeedChk(ierr);
182 ierr = CeedVectorSetValue(elem_dof, 0.0); PCeedChk(ierr);
183 CeedElemRestrictionApply(rstrin, CEED_NOTRANSPOSE, index_vec,
184 elem_dof, CEED_REQUEST_IMMEDIATE); PCeedChk(ierr);
185 const CeedScalar * elem_dof_a;
186 ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a);
187 PCeedChk(ierr);
188 ierr = CeedVectorDestroy(&index_vec); PCeedChk(ierr);
189
190 // loop over elements and put in SparseMatrix
191 // SparseMatrix * out = new SparseMatrix(nnodes, nnodes);
192 MFEM_ASSERT(out->Height() == nnodes, "Sizes don't match!");
193 MFEM_ASSERT(out->Width() == nnodes, "Sizes don't match!");
194 const CeedScalar *interpin, *gradin;
195 ierr = CeedBasisGetInterp(basisin, &interpin); PCeedChk(ierr);
196 ierr = CeedBasisGetGrad(basisin, &gradin); PCeedChk(ierr);
197
198 const CeedScalar * assembledqfarray;
199 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray);
200 PCeedChk(ierr);
201
202 CeedInt layout[3];
203 ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); PCeedChk(ierr);
204 ierr = CeedElemRestrictionDestroy(&rstr_q); PCeedChk(ierr);
205
206 // enforce structurally symmetric for later elimination
207 const int skip_zeros = 0;
208 MFEM_ASSERT(numemodein == numemodeout,
209 "Ceed full assembly not implemented for this case.");
210 for (int e = 0; e < nelem; ++e)
211 {
212 // get Array<int> for use in SparseMatrix::AddSubMatrix()
213 Array<int> rows(elemsize);
214 for (int i = 0; i < elemsize; ++i)
215 {
216 rows[i] = elem_dof_a[e * elemsize + i];
217 }
218
219 // form element matrix itself
220 DenseMatrix Bmat(nqpts * numemodein, elemsize);
221 Bmat = 0.0;
222 // Store block-diagonal D matrix as collection of small dense blocks
223 DenseTensor Dmat(numemodeout, numemodein, nqpts);
224 Dmat = 0.0;
225 DenseMatrix elem_mat(elemsize, elemsize);
226 elem_mat = 0.0;
227 for (int q = 0; q < nqpts; ++q)
228 {
229 for (int n = 0; n < elemsize; ++n)
230 {
231 CeedInt din = -1;
232 for (int ein = 0; ein < numemodein; ++ein)
233 {
234 if (emodein[ein] == CEED_EVAL_INTERP)
235 {
236 Bmat(numemodein * q + ein, n) += interpin[q * elemsize + n];
237 }
238 else if (emodein[ein] == CEED_EVAL_GRAD)
239 {
240 din += 1;
241 Bmat(numemodein * q + ein, n) += gradin[(din*nqpts+q) * elemsize + n];
242 }
243 else
244 {
245 MFEM_ASSERT(false, "Not implemented!");
246 }
247 }
248 }
249 for (int ei = 0; ei < numemodein; ++ei)
250 {
251 for (int ej = 0; ej < numemodein; ++ej)
252 {
253 const int comp = ei * numemodein + ej;
254 const int index = q*layout[0] + comp*layout[1] + e*layout[2];
255 Dmat(ei, ej, q) += assembledqfarray[index];
256 }
257 }
258 }
259 DenseMatrix BTD(elemsize, nqpts*numemodein);
260 // Compute B^T*D
261 BTD = 0.0;
262 for (int j=0; j<elemsize; ++j)
263 {
264 for (int q=0; q<nqpts; ++q)
265 {
266 int qq = numemodein*q;
267 for (int ei = 0; ei < numemodein; ++ei)
268 {
269 for (int ej = 0; ej < numemodein; ++ej)
270 {
271 BTD(j,qq+ei) += Bmat(qq+ej,j)*Dmat(ej,ei,q);
272 }
273 }
274 }
275 }
276
277 Mult(BTD, Bmat, elem_mat);
278
279 // put element matrix in sparsemat
280 out->AddSubMatrix(rows, rows, elem_mat, skip_zeros);
281 }
282
283 ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); PCeedChk(ierr);
284 ierr = CeedVectorDestroy(&elem_dof); PCeedChk(ierr);
285 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray);
286 PCeedChk(ierr);
287 ierr = CeedVectorDestroy(&assembledqf); PCeedChk(ierr);
288 ierr = CeedHackFree(&emodein); PCeedChk(ierr);
289 ierr = CeedHackFree(&emodeout); PCeedChk(ierr);
290
291 return 0;
292}
293
294int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
295{
296 int ierr;
297
298 CeedSize in_len, out_len;
299 ierr = CeedOperatorGetActiveVectorLengths(op, &in_len, &out_len);
300 PCeedChk(ierr);
301 const int nnodes = in_len;
302 MFEM_VERIFY(in_len == out_len, "not a square CeedOperator");
303 MFEM_VERIFY(in_len == nnodes, "size overflow");
304
305 SparseMatrix *out = new SparseMatrix(nnodes, nnodes);
306
307 bool isComposite;
308 ierr = CeedOperatorIsComposite(op, &isComposite); PCeedChk(ierr);
309 if (isComposite)
310 {
311 CeedInt numsub;
312 CeedOperator *subops;
313#if CEED_VERSION_GE(0, 10, 2)
314 CeedCompositeOperatorGetNumSub(op, &numsub);
315 ierr = CeedCompositeOperatorGetSubList(op, &subops); PCeedChk(ierr);
316#else
317 CeedOperatorGetNumSub(op, &numsub);
318 ierr = CeedOperatorGetSubList(op, &subops); PCeedChk(ierr);
319#endif
320 for (int i = 0; i < numsub; ++i)
321 {
322 ierr = CeedSingleOperatorFullAssemble(subops[i], out); PCeedChk(ierr);
323 }
324 }
325 else
326 {
327 ierr = CeedSingleOperatorFullAssemble(op, out); PCeedChk(ierr);
328 }
329 // enforce structurally symmetric for later elimination
330 const int skip_zeros = 0;
331 out->Finalize(skip_zeros);
332 *mat = out;
333
334 return 0;
335}
336
337} // namespace ceed
338
339} // namespace mfem
340
341#endif // MFEM_USE_CEED
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition densemat.hpp:1104
mfem::SparseMatrix
Data type sparse matrix.
Definition sparsemat.hpp:51
dim
int dim
Definition ex24.cpp:53
index
int index(int i, int j, int nx, int ny)
Definition life.cpp:235
mfem::ceed::CeedSingleOperatorFullAssemble
int CeedSingleOperatorFullAssemble(CeedOperator op, SparseMatrix *out)
Definition full-assembly.cpp:44
mfem::ceed::CeedHackFree
int CeedHackFree(void *p)
Definition full-assembly.cpp:37
mfem::ceed::CeedOperatorFullAssemble
int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
Assembles a CeedOperator as an mfem::SparseMatrix.
Definition full-assembly.cpp:294
mfem::ceed::CeedHackReallocArray
int CeedHackReallocArray(size_t n, size_t unit, void *p)
Definition full-assembly.cpp:26
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
mfem::out
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
doxygen 1.11.0