1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Diffusion operator example using MFEM 19 #include <ceed.h> 20 #include <mfem.hpp> 21 22 #include "bp3.h" 23 24 /// Wrapper for a diffusion CeedOperator as an mfem::Operator 25 class CeedDiffusionOperator : public mfem::Operator { 26 protected: 27 const mfem::FiniteElementSpace *fes; 28 CeedOperator build_oper, oper; 29 CeedBasis basis, mesh_basis; 30 CeedElemRestriction restr, mesh_restr, restr_i, mesh_restr_i; 31 CeedQFunction apply_qfunc, build_qfunc; 32 CeedQFunctionContext build_ctx; 33 CeedVector node_coords, qdata; 34 35 BuildContext build_ctx_data; 36 37 CeedVector u, v; 38 39 static void FESpace2Ceed(const mfem::FiniteElementSpace *fes, 40 const mfem::IntegrationRule &ir, 41 Ceed ceed, CeedBasis *basis, 42 CeedElemRestriction *restr) { 43 mfem::Mesh *mesh = fes->GetMesh(); 44 const mfem::FiniteElement *fe = fes->GetFE(0); 45 const int order = fes->GetOrder(0); 46 mfem::Array<int> dof_map; 47 switch (mesh->Dimension()) { 48 case 1: { 49 const mfem::H1_SegmentElement *h1_fe = 50 dynamic_cast<const mfem::H1_SegmentElement *>(fe); 51 MFEM_VERIFY(h1_fe, "invalid FE"); 52 h1_fe->GetDofMap().Copy(dof_map); 53 break; 54 } 55 case 2: { 56 const mfem::H1_QuadrilateralElement *h1_fe = 57 dynamic_cast<const mfem::H1_QuadrilateralElement *>(fe); 58 MFEM_VERIFY(h1_fe, "invalid FE"); 59 h1_fe->GetDofMap().Copy(dof_map); 60 break; 61 } 62 case 3: { 63 const mfem::H1_HexahedronElement *h1_fe = 64 dynamic_cast<const mfem::H1_HexahedronElement *>(fe); 65 MFEM_VERIFY(h1_fe, "invalid FE"); 66 h1_fe->GetDofMap().Copy(dof_map); 67 break; 68 } 69 } 70 const mfem::FiniteElement *fe1d = 71 fes->FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT); 72 mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints()); 73 mfem::DenseMatrix grad1d(fe1d->GetDof(), ir.GetNPoints()); 74 mfem::Vector qref1d(ir.GetNPoints()), qweight1d(ir.GetNPoints()); 75 mfem::Vector shape_i(shape1d.Height()); 76 mfem::DenseMatrix grad_i(grad1d.Height(), 1); 77 const mfem::H1_SegmentElement *h1_fe1d = 78 dynamic_cast<const mfem::H1_SegmentElement *>(fe1d); 79 MFEM_VERIFY(h1_fe1d, "invalid FE"); 80 const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap(); 81 for (int i = 0; i < ir.GetNPoints(); i++) { 82 const mfem::IntegrationPoint &ip = ir.IntPoint(i); 83 qref1d(i) = ip.x; 84 qweight1d(i) = ip.weight; 85 fe1d->CalcShape(ip, shape_i); 86 fe1d->CalcDShape(ip, grad_i); 87 for (int j = 0; j < shape1d.Height(); j++) { 88 shape1d(j,i) = shape_i(dof_map_1d[j]); 89 grad1d(j,i) = grad_i(dof_map_1d[j],0); 90 } 91 } 92 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes->GetVDim(), order+1, 93 ir.GetNPoints(), shape1d.GetData(), 94 grad1d.GetData(), qref1d.GetData(), 95 qweight1d.GetData(), basis); 96 97 const mfem::Table &el_dof = fes->GetElementToDofTable(); 98 mfem::Array<int> tp_el_dof(el_dof.Size_of_connections()); 99 for (int i = 0; i < mesh->GetNE(); i++) { 100 const int el_offset = fe->GetDof()*i; 101 for (int j = 0; j < fe->GetDof(); j++) { 102 tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset]; 103 } 104 } 105 CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(), 106 fes->GetVDim(), fes->GetNDofs(), 107 (fes->GetVDim())*(fes->GetNDofs()), 108 CEED_MEM_HOST, CEED_COPY_VALUES, 109 tp_el_dof.GetData(), restr); 110 } 111 112 public: 113 /// Constructor. Assumes @a fes is a scalar FE space. 114 CeedDiffusionOperator(Ceed ceed, const mfem::FiniteElementSpace *fes) 115 : Operator(fes->GetNDofs()), 116 fes(fes) { 117 mfem::Mesh *mesh = fes->GetMesh(); 118 const int order = fes->GetOrder(0); 119 const int ir_order = 2*(order + 2) - 1; // <----- 120 const mfem::IntegrationRule &ir = 121 mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order); 122 CeedInt nelem = mesh->GetNE(), dim = mesh->SpaceDimension(), 123 ncompx = dim, nqpts; 124 125 FESpace2Ceed(fes, ir, ceed, &basis, &restr); 126 127 const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace(); 128 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space"); 129 FESpace2Ceed(mesh_fes, ir, ceed, &mesh_basis, &mesh_restr); 130 CeedBasisGetNumQuadraturePoints(basis, &nqpts); 131 132 CeedInt strides[3] = {1, nqpts, nqpts *dim *(dim+1)/2}; 133 CeedElemRestrictionCreateStrided(ceed, nelem, nqpts, dim*(dim+1)/2, 134 dim*(dim+1)/2*nqpts*nelem, strides, 135 &restr_i); 136 137 CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &node_coords); 138 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, 139 mesh->GetNodes()->GetData()); 140 141 CeedVectorCreate(ceed, nelem*nqpts*dim*(dim+1)/2, &qdata); 142 143 // Context data to be passed to the 'f_build_diff' Q-function. 144 build_ctx_data.dim = mesh->Dimension(); 145 build_ctx_data.space_dim = dim; 146 CeedQFunctionContextCreate(ceed, &build_ctx); 147 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST, CEED_USE_POINTER, 148 sizeof(build_ctx_data), &build_ctx_data); 149 150 // Create the Q-function that builds the diff operator (i.e. computes its 151 // quadrature data) and set its context data. 152 CeedQFunctionCreateInterior(ceed, 1, f_build_diff, 153 f_build_diff_loc, &build_qfunc); 154 CeedQFunctionAddInput(build_qfunc, "dx", ncompx*dim, CEED_EVAL_GRAD); 155 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT); 156 CeedQFunctionAddOutput(build_qfunc, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 157 CeedQFunctionSetContext(build_qfunc, build_ctx); 158 159 // Create the operator that builds the quadrature data for the diff operator. 160 CeedOperatorCreate(ceed, build_qfunc, CEED_QFUNCTION_NONE, 161 CEED_QFUNCTION_NONE, &build_oper); 162 CeedOperatorSetField(build_oper, "dx", mesh_restr, mesh_basis, 163 CEED_VECTOR_ACTIVE); 164 CeedOperatorSetField(build_oper, "weights", CEED_ELEMRESTRICTION_NONE, 165 mesh_basis, CEED_VECTOR_NONE); 166 CeedOperatorSetField(build_oper, "qdata", restr_i, CEED_BASIS_COLLOCATED, 167 CEED_VECTOR_ACTIVE); 168 169 // Compute the quadrature data for the diff operator. 170 CeedOperatorApply(build_oper, node_coords, qdata, 171 CEED_REQUEST_IMMEDIATE); 172 173 // Create the Q-function that defines the action of the diff operator. 174 CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, 175 f_apply_diff_loc, &apply_qfunc); 176 CeedQFunctionAddInput(apply_qfunc, "u", dim, CEED_EVAL_GRAD); 177 CeedQFunctionAddInput(apply_qfunc, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 178 CeedQFunctionAddOutput(apply_qfunc, "v", dim, CEED_EVAL_GRAD); 179 CeedQFunctionSetContext(apply_qfunc, build_ctx); 180 181 // Create the diff operator. 182 CeedOperatorCreate(ceed, apply_qfunc, CEED_QFUNCTION_NONE, 183 CEED_QFUNCTION_NONE, &oper); 184 CeedOperatorSetField(oper, "u", restr, basis, CEED_VECTOR_ACTIVE); 185 CeedOperatorSetField(oper, "qdata", restr_i, CEED_BASIS_COLLOCATED, qdata); 186 CeedOperatorSetField(oper, "v", restr, basis, CEED_VECTOR_ACTIVE); 187 188 CeedVectorCreate(ceed, fes->GetNDofs(), &u); 189 CeedVectorCreate(ceed, fes->GetNDofs(), &v); 190 } 191 192 /// Destructor 193 ~CeedDiffusionOperator() { 194 CeedVectorDestroy(&u); 195 CeedVectorDestroy(&v); 196 CeedVectorDestroy(&qdata); 197 CeedVectorDestroy(&node_coords); 198 CeedElemRestrictionDestroy(&restr); 199 CeedElemRestrictionDestroy(&mesh_restr); 200 CeedElemRestrictionDestroy(&restr_i); 201 CeedBasisDestroy(&basis); 202 CeedBasisDestroy(&mesh_basis); 203 CeedQFunctionDestroy(&build_qfunc); 204 CeedQFunctionContextDestroy(&build_ctx); 205 CeedOperatorDestroy(&build_oper); 206 CeedQFunctionDestroy(&apply_qfunc); 207 CeedOperatorDestroy(&oper); 208 } 209 210 /// Operator action 211 virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const { 212 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_USE_POINTER, x.GetData()); 213 CeedVectorSetArray(v, CEED_MEM_HOST, CEED_USE_POINTER, y.GetData()); 214 215 CeedOperatorApply(oper, u, v, CEED_REQUEST_IMMEDIATE); 216 CeedVectorSyncArray(v, CEED_MEM_HOST); 217 } 218 }; 219