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