MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
full-assembly.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 CeedBasis basis;
85 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis); PCeedChk(ierr);
86 if (!basisin)
87 {
88 ierr = CeedBasisReferenceCopy(basis, &basisin); PCeedChk(ierr);
89 }
90#if CEED_VERSION_GE(0, 13, 0)
91 ierr = CeedBasisDestroy(&basis); PCeedChk(ierr);
92#endif
93 ierr = CeedBasisGetNumComponents(basisin, &ncomp); PCeedChk(ierr);
94 ierr = CeedBasisGetDimension(basisin, &dim); PCeedChk(ierr);
95 CeedElemRestriction rstr;
96 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr);
97 PCeedChk(ierr);
98 if (!rstrin)
99 {
100 ierr = CeedElemRestrictionReferenceCopy(rstr, &rstrin); PCeedChk(ierr);
101 }
102#if CEED_VERSION_GE(0, 13, 0)
103 ierr = CeedElemRestrictionDestroy(&rstr); PCeedChk(ierr);
104#endif
105 CeedEvalMode emode;
106 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
107 PCeedChk(ierr);
108 switch (emode)
109 {
110 case CEED_EVAL_NONE:
111 case CEED_EVAL_INTERP:
112 ierr = CeedHackRealloc(numemodein + 1, &emodein); PCeedChk(ierr);
113 emodein[numemodein] = emode;
114 numemodein += 1;
115 break;
116 case CEED_EVAL_GRAD:
117 ierr = CeedHackRealloc(numemodein + dim, &emodein); PCeedChk(ierr);
118 for (CeedInt d=0; d<dim; d++)
119 {
120 emodein[numemodein+d] = emode;
121 }
122 numemodein += dim;
123 break;
124 case CEED_EVAL_WEIGHT:
125 case CEED_EVAL_DIV:
126 case CEED_EVAL_CURL:
127 break; // Caught by QF Assembly
128 }
129 }
130#if CEED_VERSION_GE(0, 13, 0)
131 ierr = CeedVectorDestroy(&vec); PCeedChk(ierr);
132#endif
133 }
134
135 // Determine active output basis
136 ierr = CeedQFunctionGetFields(qf, &numinputfields, NULL, &numoutputfields,
137 &qffields); PCeedChk(ierr);
138 CeedInt numemodeout = 0;
139 CeedEvalMode *emodeout = NULL;
140 CeedBasis basisout = NULL;
141 CeedElemRestriction rstrout = NULL;
142 for (CeedInt i=0; i<numoutputfields; i++)
143 {
144 CeedVector vec;
145 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); PCeedChk(ierr);
146 if (vec == CEED_VECTOR_ACTIVE)
147 {
148 CeedBasis basis;
149 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis); PCeedChk(ierr);
150 if (!basisout)
151 {
152 ierr = CeedBasisReferenceCopy(basis, &basisout); PCeedChk(ierr);
153 }
154#if CEED_VERSION_GE(0, 13, 0)
155 ierr = CeedBasisDestroy(&basis); PCeedChk(ierr);
156#endif
157 CeedElemRestriction rstr;
158 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr);
159 PCeedChk(ierr);
160 if (!rstrout)
161 {
162 ierr = CeedElemRestrictionReferenceCopy(rstr, &rstrout); PCeedChk(ierr);
163 }
164#if CEED_VERSION_GE(0, 13, 0)
165 ierr = CeedElemRestrictionDestroy(&rstr); PCeedChk(ierr);
166#endif
167 CeedEvalMode emode;
168 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
169 PCeedChk(ierr);
170 switch (emode)
171 {
172 case CEED_EVAL_NONE:
173 case CEED_EVAL_INTERP:
174 ierr = CeedHackRealloc(numemodeout + 1, &emodeout); PCeedChk(ierr);
175 emodeout[numemodeout] = emode;
176 numemodeout += 1;
177 break;
178 case CEED_EVAL_GRAD:
179 ierr = CeedHackRealloc(numemodeout + dim, &emodeout); PCeedChk(ierr);
180 for (CeedInt d=0; d<dim; d++)
181 {
182 emodeout[numemodeout+d] = emode;
183 }
184 numemodeout += dim;
185 break;
186 case CEED_EVAL_WEIGHT:
187 case CEED_EVAL_DIV:
188 case CEED_EVAL_CURL:
189 break; // Caught by QF Assembly
190 }
191 }
192#if CEED_VERSION_GE(0, 13, 0)
193 ierr = CeedVectorDestroy(&vec); PCeedChk(ierr);
194#endif
195 }
196
197 CeedInt nelem, elemsize, nqpts;
198 CeedSize nnodes;
199 ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); PCeedChk(ierr);
200 ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); PCeedChk(ierr);
201 ierr = CeedElemRestrictionGetLVectorSize(rstrin, &nnodes); PCeedChk(ierr);
202 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); PCeedChk(ierr);
203
204 // Determine elem_dof relation
205 CeedVector index_vec;
206 ierr = CeedVectorCreate(ceed, nnodes, &index_vec); PCeedChk(ierr);
207 CeedScalar *array;
208 ierr = CeedVectorGetArrayWrite(index_vec, CEED_MEM_HOST, &array);
209 PCeedChk(ierr);
210 for (CeedSize i = 0; i < nnodes; ++i)
211 {
212 array[i] = i;
213 }
214 ierr = CeedVectorRestoreArray(index_vec, &array); PCeedChk(ierr);
215 CeedVector elem_dof;
216 ierr = CeedVectorCreate(ceed, nelem * elemsize, &elem_dof); PCeedChk(ierr);
217 ierr = CeedVectorSetValue(elem_dof, 0.0); PCeedChk(ierr);
218 CeedElemRestrictionApply(rstrin, CEED_NOTRANSPOSE, index_vec,
219 elem_dof, CEED_REQUEST_IMMEDIATE); PCeedChk(ierr);
220 const CeedScalar * elem_dof_a;
221 ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a);
222 PCeedChk(ierr);
223 ierr = CeedVectorDestroy(&index_vec); PCeedChk(ierr);
224
225 // loop over elements and put in SparseMatrix
226 // SparseMatrix * out = new SparseMatrix(nnodes, nnodes);
227 MFEM_ASSERT(out->Height() == nnodes, "Sizes don't match!");
228 MFEM_ASSERT(out->Width() == nnodes, "Sizes don't match!");
229 const CeedScalar *interpin, *gradin;
230 ierr = CeedBasisGetInterp(basisin, &interpin); PCeedChk(ierr);
231 ierr = CeedBasisGetGrad(basisin, &gradin); PCeedChk(ierr);
232
233 const CeedScalar * assembledqfarray;
234 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray);
235 PCeedChk(ierr);
236
237 CeedInt layout[3];
238#if CEED_VERSION_GE(0, 13, 0)
239 ierr = CeedElemRestrictionGetELayout(rstr_q, layout); PCeedChk(ierr);
240#else
241 ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); PCeedChk(ierr);
242#endif
243 ierr = CeedElemRestrictionDestroy(&rstr_q); PCeedChk(ierr);
244
245 // enforce structurally symmetric for later elimination
246 const int skip_zeros = 0;
247 MFEM_ASSERT(numemodein == numemodeout,
248 "Ceed full assembly not implemented for this case.");
249 for (int e = 0; e < nelem; ++e)
250 {
251 // get Array<int> for use in SparseMatrix::AddSubMatrix()
252 Array<int> rows(elemsize);
253 for (int i = 0; i < elemsize; ++i)
254 {
255 rows[i] = elem_dof_a[e * elemsize + i];
256 }
257
258 // form element matrix itself
259 DenseMatrix Bmat(nqpts * numemodein, elemsize);
260 Bmat = 0.0;
261 // Store block-diagonal D matrix as collection of small dense blocks
262 DenseTensor Dmat(numemodeout, numemodein, nqpts);
263 Dmat = 0.0;
264 DenseMatrix elem_mat(elemsize, elemsize);
265 elem_mat = 0.0;
266 for (int q = 0; q < nqpts; ++q)
267 {
268 for (int n = 0; n < elemsize; ++n)
269 {
270 CeedInt din = -1;
271 for (int ein = 0; ein < numemodein; ++ein)
272 {
273 if (emodein[ein] == CEED_EVAL_INTERP)
274 {
275 Bmat(numemodein * q + ein, n) += interpin[q * elemsize + n];
276 }
277 else if (emodein[ein] == CEED_EVAL_GRAD)
278 {
279 din += 1;
280 Bmat(numemodein * q + ein, n) += gradin[(din*nqpts+q) * elemsize + n];
281 }
282 else
283 {
284 MFEM_ASSERT(false, "Not implemented!");
285 }
286 }
287 }
288 for (int ei = 0; ei < numemodein; ++ei)
289 {
290 for (int ej = 0; ej < numemodein; ++ej)
291 {
292 const int comp = ei * numemodein + ej;
293 const int index = q*layout[0] + comp*layout[1] + e*layout[2];
294 Dmat(ei, ej, q) += assembledqfarray[index];
295 }
296 }
297 }
298 DenseMatrix BTD(elemsize, nqpts*numemodein);
299 // Compute B^T*D
300 BTD = 0.0;
301 for (int j=0; j<elemsize; ++j)
302 {
303 for (int q=0; q<nqpts; ++q)
304 {
305 int qq = numemodein*q;
306 for (int ei = 0; ei < numemodein; ++ei)
307 {
308 for (int ej = 0; ej < numemodein; ++ej)
309 {
310 BTD(j,qq+ei) += Bmat(qq+ej,j)*Dmat(ej,ei,q);
311 }
312 }
313 }
314 }
315
316 Mult(BTD, Bmat, elem_mat);
317
318 // put element matrix in sparsemat
319 out->AddSubMatrix(rows, rows, elem_mat, skip_zeros);
320 }
321
322 ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); PCeedChk(ierr);
323 ierr = CeedVectorDestroy(&elem_dof); PCeedChk(ierr);
324 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray);
325 PCeedChk(ierr);
326 ierr = CeedVectorDestroy(&assembledqf); PCeedChk(ierr);
327 ierr = CeedElemRestrictionDestroy(&rstrin); PCeedChk(ierr);
328 ierr = CeedElemRestrictionDestroy(&rstrout); PCeedChk(ierr);
329 ierr = CeedBasisDestroy(&basisin); PCeedChk(ierr);
330 ierr = CeedBasisDestroy(&basisout); PCeedChk(ierr);
331 ierr = CeedHackFree(&emodein); PCeedChk(ierr);
332 ierr = CeedHackFree(&emodeout); PCeedChk(ierr);
333
334 return 0;
335}
336
337int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
338{
339 int ierr;
340
341 CeedSize in_len, out_len;
342 ierr = CeedOperatorGetActiveVectorLengths(op, &in_len, &out_len);
343 PCeedChk(ierr);
344 const int nnodes = in_len;
345 MFEM_VERIFY(in_len == out_len, "not a square CeedOperator");
346 MFEM_VERIFY(in_len == nnodes, "size overflow");
347
348 SparseMatrix *out = new SparseMatrix(nnodes, nnodes);
349
350 bool isComposite;
351 ierr = CeedOperatorIsComposite(op, &isComposite); PCeedChk(ierr);
352 if (isComposite)
353 {
354 CeedInt numsub;
355 CeedOperator *subops;
356 ierr = CeedOperatorCompositeGetNumSub(op, &numsub); PCeedChk(ierr);
357 ierr = CeedOperatorCompositeGetSubList(op, &subops); PCeedChk(ierr);
358 for (int i = 0; i < numsub; ++i)
359 {
360 ierr = CeedSingleOperatorFullAssemble(subops[i], out); PCeedChk(ierr);
361 }
362 }
363 else
364 {
365 ierr = CeedSingleOperatorFullAssemble(op, out); PCeedChk(ierr);
366 }
367 // enforce structurally symmetric for later elimination
368 const int skip_zeros = 0;
369 out->Finalize(skip_zeros);
370 *mat = out;
371
372 return 0;
373}
374
375} // namespace ceed
376
377} // namespace mfem
378
379#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:1111
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:236
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:337
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:505
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
Generated on Thu Dec 11 2025 11:52:53 for MFEM by doxygen 1.11.0