1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Diffusion operator example using MFEM 10 11 #include <ceed.h> 12 13 #include <mfem.hpp> 14 15 #include "bp3.h" 16 17 /// Wrapper for a diffusion CeedOperator as an mfem::Operator 18 class CeedDiffusionOperator : public mfem::Operator { 19 protected: 20 const mfem::FiniteElementSpace *fes; 21 CeedOperator build_oper, oper; 22 CeedBasis basis, mesh_basis; 23 CeedElemRestriction restr, mesh_restr, restr_i, mesh_restr_i; 24 CeedQFunction apply_qfunc, build_qfunc; 25 CeedQFunctionContext build_ctx; 26 CeedVector node_coords, qdata; 27 28 BuildContext build_ctx_data; 29 30 CeedVector u, v; 31 32 static void FESpace2Ceed(const mfem::FiniteElementSpace *fes, const mfem::IntegrationRule &ir, Ceed ceed, CeedBasis *basis, 33 CeedElemRestriction *restr) { 34 mfem::Mesh *mesh = fes->GetMesh(); 35 const mfem::FiniteElement *fe = fes->GetFE(0); 36 const int order = fes->GetOrder(0); 37 mfem::Array<int> dof_map; 38 switch (mesh->Dimension()) { 39 case 1: { 40 const mfem::H1_SegmentElement *h1_fe = dynamic_cast<const mfem::H1_SegmentElement *>(fe); 41 MFEM_VERIFY(h1_fe, "invalid FE"); 42 h1_fe->GetDofMap().Copy(dof_map); 43 break; 44 } 45 case 2: { 46 const mfem::H1_QuadrilateralElement *h1_fe = dynamic_cast<const mfem::H1_QuadrilateralElement *>(fe); 47 MFEM_VERIFY(h1_fe, "invalid FE"); 48 h1_fe->GetDofMap().Copy(dof_map); 49 break; 50 } 51 case 3: { 52 const mfem::H1_HexahedronElement *h1_fe = dynamic_cast<const mfem::H1_HexahedronElement *>(fe); 53 MFEM_VERIFY(h1_fe, "invalid FE"); 54 h1_fe->GetDofMap().Copy(dof_map); 55 break; 56 } 57 } 58 const mfem::FiniteElement *fe1d = fes->FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT); 59 mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints()); 60 mfem::DenseMatrix grad_1d(fe1d->GetDof(), ir.GetNPoints()); 61 mfem::Vector q_ref_1d(ir.GetNPoints()), q_weight_1d(ir.GetNPoints()); 62 mfem::Vector shape_i(shape1d.Height()); 63 mfem::DenseMatrix grad_i(grad_1d.Height(), 1); 64 const mfem::H1_SegmentElement *h1_fe1d = dynamic_cast<const mfem::H1_SegmentElement *>(fe1d); 65 MFEM_VERIFY(h1_fe1d, "invalid FE"); 66 const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap(); 67 for (int i = 0; i < ir.GetNPoints(); i++) { 68 const mfem::IntegrationPoint &ip = ir.IntPoint(i); 69 q_ref_1d(i) = ip.x; 70 q_weight_1d(i) = ip.weight; 71 fe1d->CalcShape(ip, shape_i); 72 fe1d->CalcDShape(ip, grad_i); 73 for (int j = 0; j < shape1d.Height(); j++) { 74 shape1d(j, i) = shape_i(dof_map_1d[j]); 75 grad_1d(j, i) = grad_i(dof_map_1d[j], 0); 76 } 77 } 78 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes->GetVDim(), order + 1, ir.GetNPoints(), shape1d.GetData(), grad_1d.GetData(), 79 q_ref_1d.GetData(), q_weight_1d.GetData(), basis); 80 81 const mfem::Table &el_dof = fes->GetElementToDofTable(); 82 mfem::Array<int> tp_el_dof(el_dof.Size_of_connections()); 83 for (int i = 0; i < mesh->GetNE(); i++) { 84 const int el_offset = fe->GetDof() * i; 85 for (int j = 0; j < fe->GetDof(); j++) { 86 tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset]; 87 } 88 } 89 CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(), fes->GetVDim(), fes->GetNDofs(), (fes->GetVDim()) * (fes->GetNDofs()), CEED_MEM_HOST, 90 CEED_COPY_VALUES, tp_el_dof.GetData(), restr); 91 } 92 93 public: 94 /// Constructor. Assumes @a fes is a scalar FE space. 95 CeedDiffusionOperator(Ceed ceed, const mfem::FiniteElementSpace *fes) : Operator(fes->GetNDofs()), fes(fes) { 96 mfem::Mesh *mesh = fes->GetMesh(); 97 const int order = fes->GetOrder(0); 98 const int ir_order = 2 * (order + 2) - 1; // <----- 99 const mfem::IntegrationRule &ir = mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order); 100 CeedInt num_elem = mesh->GetNE(), dim = mesh->SpaceDimension(), ncompx = dim, nqpts; 101 102 FESpace2Ceed(fes, ir, ceed, &basis, &restr); 103 104 const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace(); 105 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space"); 106 FESpace2Ceed(mesh_fes, ir, ceed, &mesh_basis, &mesh_restr); 107 CeedBasisGetNumQuadraturePoints(basis, &nqpts); 108 109 CeedInt strides[3] = {1, nqpts, nqpts * dim * (dim + 1) / 2}; 110 CeedElemRestrictionCreateStrided(ceed, num_elem, nqpts, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * nqpts * num_elem, strides, &restr_i); 111 112 CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &node_coords); 113 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, mesh->GetNodes()->GetData()); 114 115 CeedVectorCreate(ceed, num_elem * nqpts * dim * (dim + 1) / 2, &qdata); 116 117 // Context data to be passed to the 'f_build_diff' Q-function. 118 build_ctx_data.dim = mesh->Dimension(); 119 build_ctx_data.space_dim = dim; 120 CeedQFunctionContextCreate(ceed, &build_ctx); 121 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(build_ctx_data), &build_ctx_data); 122 123 // Create the Q-function that builds the diff operator (i.e. computes its quadrature data) and set its context data. 124 CeedQFunctionCreateInterior(ceed, 1, f_build_diff, f_build_diff_loc, &build_qfunc); 125 CeedQFunctionAddInput(build_qfunc, "dx", ncompx * dim, CEED_EVAL_GRAD); 126 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT); 127 CeedQFunctionAddOutput(build_qfunc, "qdata", dim * (dim + 1) / 2, CEED_EVAL_NONE); 128 CeedQFunctionSetContext(build_qfunc, build_ctx); 129 130 // Create the operator that builds the quadrature data for the diff operator. 131 CeedOperatorCreate(ceed, build_qfunc, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &build_oper); 132 CeedOperatorSetField(build_oper, "dx", mesh_restr, mesh_basis, CEED_VECTOR_ACTIVE); 133 CeedOperatorSetField(build_oper, "weights", CEED_ELEMRESTRICTION_NONE, mesh_basis, CEED_VECTOR_NONE); 134 CeedOperatorSetField(build_oper, "qdata", restr_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 135 136 // Compute the quadrature data for the diff operator. 137 CeedOperatorApply(build_oper, node_coords, qdata, CEED_REQUEST_IMMEDIATE); 138 139 // Create the Q-function that defines the action of the diff operator. 140 CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, f_apply_diff_loc, &apply_qfunc); 141 CeedQFunctionAddInput(apply_qfunc, "u", dim, CEED_EVAL_GRAD); 142 CeedQFunctionAddInput(apply_qfunc, "qdata", dim * (dim + 1) / 2, CEED_EVAL_NONE); 143 CeedQFunctionAddOutput(apply_qfunc, "v", dim, CEED_EVAL_GRAD); 144 CeedQFunctionSetContext(apply_qfunc, build_ctx); 145 146 // Create the diff operator. 147 CeedOperatorCreate(ceed, apply_qfunc, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &oper); 148 CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE); 149 CeedOperatorSetField(oper, "qdata", restr_i, CEED_BASIS_NONE, qdata); 150 CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE); 151 152 CeedVectorCreate(ceed, fes->GetNDofs(), &u); 153 CeedVectorCreate(ceed, fes->GetNDofs(), &v); 154 } 155 156 /// Destructor 157 ~CeedDiffusionOperator() { 158 CeedVectorDestroy(&u); 159 CeedVectorDestroy(&v); 160 CeedVectorDestroy(&qdata); 161 CeedVectorDestroy(&node_coords); 162 CeedElemRestrictionDestroy(&restr); 163 CeedElemRestrictionDestroy(&mesh_restr); 164 CeedElemRestrictionDestroy(&restr_i); 165 CeedBasisDestroy(&basis); 166 CeedBasisDestroy(&mesh_basis); 167 CeedQFunctionDestroy(&build_qfunc); 168 CeedQFunctionContextDestroy(&build_ctx); 169 CeedOperatorDestroy(&build_oper); 170 CeedQFunctionDestroy(&apply_qfunc); 171 CeedOperatorDestroy(&oper); 172 } 173 174 /// Operator action 175 virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const { 176 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_USE_POINTER, x.GetData()); 177 CeedVectorSetArray(v, CEED_MEM_HOST, CEED_USE_POINTER, y.GetData()); 178 179 CeedOperatorApply(oper, u, v, CEED_REQUEST_IMMEDIATE); 180 CeedVectorSyncArray(v, CEED_MEM_HOST); 181 } 182 }; 183