// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. #include "full-assembly.hpp" #include "../../../linalg/sparsemat.hpp" #include "../interface/util.hpp" #include "../interface/ceed.hpp" #ifdef MFEM_USE_CEED namespace mfem { namespace ceed { int CeedHackReallocArray(size_t n, size_t unit, void *p) { *(void **)p = realloc(*(void **)p, n*unit); if (n && unit && !*(void **)p) return CeedError(NULL, 1, "realloc failed to allocate %zd members of size " "%zd\n", n, unit); return 0; } #define CeedHackRealloc(n, p) CeedHackReallocArray((n), sizeof(**(p)), p) int CeedHackFree(void *p) { free(*(void **)p); *(void **)p = NULL; return 0; } int CeedSingleOperatorFullAssemble(CeedOperator op, SparseMatrix *out) { int ierr; Ceed ceed; ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); // Assemble QFunction CeedQFunction qf; ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); CeedInt numinputfields, numoutputfields; CeedChk(ierr); CeedVector assembledqf; CeedElemRestriction rstr_q; ierr = CeedOperatorLinearAssembleQFunction( op, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); CeedSize qflength; ierr = CeedVectorGetLength(assembledqf, &qflength); CeedChk(ierr); CeedOperatorField *input_fields; CeedOperatorField *output_fields; ierr = CeedOperatorGetFields(op, &numinputfields, &input_fields, &numoutputfields, &output_fields); CeedChk(ierr); // Determine active input basis CeedQFunctionField *qffields; ierr = CeedQFunctionGetFields(qf, &numinputfields, &qffields, &numoutputfields, NULL); CeedChk(ierr); CeedInt numemodein = 0, ncomp, dim = 1; CeedEvalMode *emodein = NULL; CeedBasis basisin = NULL; CeedElemRestriction rstrin = NULL; for (CeedInt i=0; iHeight() == nnodes, "Sizes don't match!"); MFEM_ASSERT(out->Width() == nnodes, "Sizes don't match!"); const CeedScalar *interpin, *gradin; ierr = CeedBasisGetInterp(basisin, &interpin); CeedChk(ierr); ierr = CeedBasisGetGrad(basisin, &gradin); CeedChk(ierr); const CeedScalar * assembledqfarray; ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray); CeedChk(ierr); CeedInt layout[3]; ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); CeedChk(ierr); ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr); // enforce structurally symmetric for later elimination const int skip_zeros = 0; MFEM_ASSERT(numemodein == numemodeout, "Ceed full assembly not implemented for this case."); for (int e = 0; e < nelem; ++e) { // get Array for use in SparseMatrix::AddSubMatrix() Array rows(elemsize); for (int i = 0; i < elemsize; ++i) { rows[i] = elem_dof_a[e * elemsize + i]; } // form element matrix itself DenseMatrix Bmat(nqpts * numemodein, elemsize); Bmat = 0.0; // Store block-diagonal D matrix as collection of small dense blocks DenseTensor Dmat(numemodeout, numemodein, nqpts); Dmat = 0.0; DenseMatrix elem_mat(elemsize, elemsize); elem_mat = 0.0; for (int q = 0; q < nqpts; ++q) { for (int n = 0; n < elemsize; ++n) { CeedInt din = -1; for (int ein = 0; ein < numemodein; ++ein) { if (emodein[ein] == CEED_EVAL_INTERP) { Bmat(numemodein * q + ein, n) += interpin[q * elemsize + n]; } else if (emodein[ein] == CEED_EVAL_GRAD) { din += 1; Bmat(numemodein * q + ein, n) += gradin[(din*nqpts+q) * elemsize + n]; } else { MFEM_ASSERT(false, "Not implemented!"); } } } for (int ei = 0; ei < numemodein; ++ei) { for (int ej = 0; ej < numemodein; ++ej) { const int comp = ei * numemodein + ej; const int index = q*layout[0] + comp*layout[1] + e*layout[2]; Dmat(ei, ej, q) += assembledqfarray[index]; } } } DenseMatrix BTD(elemsize, nqpts*numemodein); // Compute B^T*D BTD = 0.0; for (int j=0; jAddSubMatrix(rows, rows, elem_mat, skip_zeros); } ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); CeedChk(ierr); ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr); ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); CeedChk(ierr); ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr); ierr = CeedHackFree(&emodein); CeedChk(ierr); ierr = CeedHackFree(&emodeout); CeedChk(ierr); return 0; } int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat) { int ierr; CeedSize in_len, out_len; ierr = CeedOperatorGetActiveVectorLengths(op, &in_len, &out_len); CeedChk(ierr); const int nnodes = in_len; MFEM_VERIFY(in_len == out_len, "not a square CeedOperator"); MFEM_VERIFY(in_len == nnodes, "size overflow"); SparseMatrix *out = new SparseMatrix(nnodes, nnodes); bool isComposite; ierr = CeedOperatorIsComposite(op, &isComposite); CeedChk(ierr); if (isComposite) { CeedInt numsub; CeedOperator *subops; #if CEED_VERSION_GE(0, 10, 2) CeedCompositeOperatorGetNumSub(op, &numsub); ierr = CeedCompositeOperatorGetSubList(op, &subops); CeedChk(ierr); #else CeedOperatorGetNumSub(op, &numsub); ierr = CeedOperatorGetSubList(op, &subops); CeedChk(ierr); #endif for (int i = 0; i < numsub; ++i) { ierr = CeedSingleOperatorFullAssemble(subops[i], out); CeedChk(ierr); } } else { ierr = CeedSingleOperatorFullAssemble(op, out); CeedChk(ierr); } // enforce structurally symmetric for later elimination const int skip_zeros = 0; out->Finalize(skip_zeros); *mat = out; return 0; } } // namespace ceed } // namespace mfem #endif