MFEM  v4.5.2
Finite element discretization library
full-assembly.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 "full-assembly.hpp"
13 
14 #include "../../../linalg/sparsemat.hpp"
15 #include "../interface/util.hpp"
16 #include "../interface/ceed.hpp"
17 
18 #ifdef MFEM_USE_CEED
19 
20 namespace mfem
21 {
22 
23 namespace ceed
24 {
25 
26 int 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 
37 int CeedHackFree(void *p)
38 {
39  free(*(void **)p);
40  *(void **)p = NULL;
41  return 0;
42 }
43 
45 {
46  int ierr;
47  Ceed ceed;
48  ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
49 
50  // Assemble QFunction
51  CeedQFunction qf;
52  ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
53  CeedInt numinputfields, numoutputfields;
54  CeedChk(ierr);
55  CeedVector assembledqf;
56  CeedElemRestriction rstr_q;
57  ierr = CeedOperatorLinearAssembleQFunction(
58  op, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
59 
60  CeedSize qflength;
61  ierr = CeedVectorGetLength(assembledqf, &qflength); CeedChk(ierr);
62 
63  CeedOperatorField *input_fields;
64  CeedOperatorField *output_fields;
65  ierr = CeedOperatorGetFields(op, &numinputfields, &input_fields,
66  &numoutputfields, &output_fields);
67  CeedChk(ierr);
68 
69  // Determine active input basis
70  CeedQFunctionField *qffields;
71  ierr = CeedQFunctionGetFields(qf, &numinputfields, &qffields,
72  &numoutputfields, NULL);
73  CeedChk(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); CeedChk(ierr);
82  if (vec == CEED_VECTOR_ACTIVE)
83  {
84  ierr = CeedOperatorFieldGetBasis(input_fields[i], &basisin);
85  CeedChk(ierr);
86  ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
87  ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
88  ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstrin);
89  CeedChk(ierr);
90  CeedEvalMode emode;
91  ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
92  CeedChk(ierr);
93  switch (emode)
94  {
95  case CEED_EVAL_NONE:
96  case CEED_EVAL_INTERP:
97  ierr = CeedHackRealloc(numemodein + 1, &emodein); CeedChk(ierr);
98  emodein[numemodein] = emode;
99  numemodein += 1;
100  break;
101  case CEED_EVAL_GRAD:
102  ierr = CeedHackRealloc(numemodein + dim, &emodein); CeedChk(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); CeedChk(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); CeedChk(ierr);
128  if (vec == CEED_VECTOR_ACTIVE)
129  {
130  ierr = CeedOperatorFieldGetBasis(output_fields[i], &basisout);
131  CeedChk(ierr);
132  ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstrout);
133  CeedChk(ierr);
134  CeedChk(ierr);
135  CeedEvalMode emode;
136  ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
137  CeedChk(ierr);
138  switch (emode)
139  {
140  case CEED_EVAL_NONE:
141  case CEED_EVAL_INTERP:
142  ierr = CeedHackRealloc(numemodeout + 1, &emodeout); CeedChk(ierr);
143  emodeout[numemodeout] = emode;
144  numemodeout += 1;
145  break;
146  case CEED_EVAL_GRAD:
147  ierr = CeedHackRealloc(numemodeout + dim, &emodeout); CeedChk(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); CeedChk(ierr);
165  ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); CeedChk(ierr);
166  ierr = CeedElemRestrictionGetLVectorSize(rstrin, &nnodes); CeedChk(ierr);
167  ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChk(ierr);
168 
169  // Determine elem_dof relation
170  CeedVector index_vec;
171  ierr = CeedVectorCreate(ceed, nnodes, &index_vec); CeedChk(ierr);
172  CeedScalar *array;
173  ierr = CeedVectorGetArrayWrite(index_vec, CEED_MEM_HOST, &array);
174  CeedChk(ierr);
175  for (CeedSize i = 0; i < nnodes; ++i)
176  {
177  array[i] = i;
178  }
179  ierr = CeedVectorRestoreArray(index_vec, &array); CeedChk(ierr);
180  CeedVector elem_dof;
181  ierr = CeedVectorCreate(ceed, nelem * elemsize, &elem_dof); CeedChk(ierr);
182  ierr = CeedVectorSetValue(elem_dof, 0.0); CeedChk(ierr);
183  CeedElemRestrictionApply(rstrin, CEED_NOTRANSPOSE, index_vec,
184  elem_dof, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
185  const CeedScalar * elem_dof_a;
186  ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a);
187  CeedChk(ierr);
188  ierr = CeedVectorDestroy(&index_vec); CeedChk(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); CeedChk(ierr);
196  ierr = CeedBasisGetGrad(basisin, &gradin); CeedChk(ierr);
197 
198  const CeedScalar * assembledqfarray;
199  ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray);
200  CeedChk(ierr);
201 
202  CeedInt layout[3];
203  ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); CeedChk(ierr);
204  ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(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); CeedChk(ierr);
284  ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr);
285  ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray);
286  CeedChk(ierr);
287  ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr);
288  ierr = CeedHackFree(&emodein); CeedChk(ierr);
289  ierr = CeedHackFree(&emodeout); CeedChk(ierr);
290 
291  return 0;
292 }
293 
294 int 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  CeedChk(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); CeedChk(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); CeedChk(ierr);
316 #else
317  CeedOperatorGetNumSub(op, &numsub);
318  ierr = CeedOperatorGetSubList(op, &subops); CeedChk(ierr);
319 #endif
320  for (int i = 0; i < numsub; ++i)
321  {
322  ierr = CeedSingleOperatorFullAssemble(subops[i], out); CeedChk(ierr);
323  }
324  }
325  else
326  {
327  ierr = CeedSingleOperatorFullAssemble(op, out); CeedChk(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
int CeedHackReallocArray(size_t n, size_t unit, void *p)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Data type sparse matrix.
Definition: sparsemat.hpp:50
int CeedSingleOperatorFullAssemble(CeedOperator op, SparseMatrix *out)
int CeedHackFree(void *p)
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
int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
Assembles a CeedOperator as an mfem::SparseMatrix.
int dim
Definition: ex24.cpp:53
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:978