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 /// MFEM mass operator based on libCEED 19 20 #include <ceed.h> 21 #include <mfem.hpp> 22 23 /// A structure used to pass additional data to f_build_mass 24 struct BuildContext { CeedInt dim, space_dim; }; 25 26 /// libCEED Q-function for building quadrature data for a mass operator 27 static int f_build_mass(void *ctx, void *qdata, CeedInt Q, 28 const CeedScalar *const *u, CeedScalar *const *v) { 29 // u[1] is Jacobians, size (Q x nc x dim) with column-major layout 30 // u[4] is quadrature weights, size (Q) 31 BuildContext *bc = (BuildContext*)ctx; 32 CeedScalar *qd = (CeedScalar*)qdata; 33 const CeedScalar *J = u[1], *qw = u[4]; 34 switch (bc->dim + 10*bc->space_dim) { 35 case 11: 36 for (CeedInt i=0; i<Q; i++) { 37 qd[i] = J[i] * qw[i]; 38 } 39 break; 40 case 22: 41 for (CeedInt i=0; i<Q; i++) { 42 // 0 2 43 // 1 3 44 qd[i] = (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * qw[i]; 45 } 46 break; 47 case 33: 48 for (CeedInt i=0; i<Q; i++) { 49 // 0 3 6 50 // 1 4 7 51 // 2 5 8 52 qd[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) - 53 J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) + 54 J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * qw[i]; 55 } 56 break; 57 default: 58 return CeedError(NULL, 1, "dim=%d, space_dim=%d is not supported", 59 bc->dim, bc->space_dim); 60 } 61 return 0; 62 } 63 64 /// libCEED Q-function for applying a mass operator 65 static int f_apply_mass(void *ctx, void *qdata, CeedInt Q, 66 const CeedScalar *const *u, CeedScalar *const *v) { 67 const CeedScalar *w = (const CeedScalar*)qdata; 68 for (CeedInt i=0; i<Q; i++) { 69 v[0][i] = w[i] * u[0][i]; 70 } 71 return 0; 72 } 73 74 /// Wrapper for a mass CeedOperator as an mfem::Operator 75 class CeedMassOperator : public mfem::Operator { 76 protected: 77 const mfem::FiniteElementSpace *fes; 78 CeedOperator build_oper, oper; 79 CeedBasis basis, mesh_basis; 80 CeedElemRestriction restr, mesh_restr; 81 CeedQFunction apply_qfunc, build_qfunc; 82 CeedVector node_coords, qdata; 83 84 BuildContext build_ctx; 85 86 CeedVector u, v; 87 88 static void FESpace2Ceed(const mfem::FiniteElementSpace *fes, 89 const mfem::IntegrationRule &ir, 90 Ceed ceed, CeedBasis *basis, 91 CeedElemRestriction *restr) { 92 mfem::Mesh *mesh = fes->GetMesh(); 93 const mfem::FiniteElement *fe = fes->GetFE(0); 94 const int order = fes->GetOrder(0); 95 mfem::Array<int> dof_map; 96 switch (mesh->Dimension()) { 97 case 1: { 98 const mfem::H1_SegmentElement *h1_fe = 99 dynamic_cast<const mfem::H1_SegmentElement*>(fe); 100 MFEM_VERIFY(h1_fe, "invalid FE"); 101 h1_fe->GetDofMap().Copy(dof_map); 102 break; 103 } 104 case 2: { 105 const mfem::H1_QuadrilateralElement *h1_fe = 106 dynamic_cast<const mfem::H1_QuadrilateralElement*>(fe); 107 MFEM_VERIFY(h1_fe, "invalid FE"); 108 h1_fe->GetDofMap().Copy(dof_map); 109 break; 110 } 111 case 3: { 112 const mfem::H1_HexahedronElement *h1_fe = 113 dynamic_cast<const mfem::H1_HexahedronElement*>(fe); 114 MFEM_VERIFY(h1_fe, "invalid FE"); 115 h1_fe->GetDofMap().Copy(dof_map); 116 break; 117 } 118 } 119 const mfem::FiniteElement *fe1d = 120 fes->FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT); 121 mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints()); 122 mfem::DenseMatrix grad1d(fe1d->GetDof(), ir.GetNPoints()); 123 mfem::Vector qref1d(ir.GetNPoints()), qweight1d(ir.GetNPoints()); 124 mfem::Vector shape_i(shape1d.Height()); 125 mfem::DenseMatrix grad_i(grad1d.Height(), 1); 126 const mfem::H1_SegmentElement *h1_fe1d = 127 dynamic_cast<const mfem::H1_SegmentElement*>(fe1d); 128 MFEM_VERIFY(h1_fe1d, "invalid FE"); 129 const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap(); 130 for (int i = 0; i < ir.GetNPoints(); i++) { 131 const mfem::IntegrationPoint &ip = ir.IntPoint(i); 132 qref1d(i) = ip.x; 133 qweight1d(i) = ip.weight; 134 fe1d->CalcShape(ip, shape_i); 135 fe1d->CalcDShape(ip, grad_i); 136 for (int j = 0; j < shape1d.Height(); j++) { 137 shape1d(j,i) = shape_i(dof_map_1d[j]); 138 grad1d(j,i) = grad_i(dof_map_1d[j],0); 139 } 140 } 141 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes->GetVDim(), order+1, 142 ir.GetNPoints(), shape1d.GetData(), 143 grad1d.GetData(), qref1d.GetData(), 144 qweight1d.GetData(), basis); 145 146 const mfem::Table &el_dof = fes->GetElementToDofTable(); 147 mfem::Array<int> tp_el_dof(el_dof.Size_of_connections()); 148 for (int i = 0; i < mesh->GetNE(); i++) { 149 const int el_offset = fe->GetDof()*i; 150 for (int j = 0; j < fe->GetDof(); j++) { 151 tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset]; 152 } 153 } 154 CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(), 155 fes->GetNDofs(), CEED_MEM_HOST, CEED_COPY_VALUES, 156 tp_el_dof.GetData(), restr); 157 } 158 159 public: 160 /// Constructor. Assumes @a fes is a scalar FE space. 161 CeedMassOperator(Ceed ceed, const mfem::FiniteElementSpace *fes) 162 : Operator(fes->GetNDofs()), 163 fes(fes) { 164 mfem::Mesh *mesh = fes->GetMesh(); 165 const int order = fes->GetOrder(0); 166 const int ir_order = 2*(order + 2) - 1; // <----- 167 const mfem::IntegrationRule &ir = 168 mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order); 169 170 FESpace2Ceed(fes, ir, ceed, &basis, &restr); 171 172 const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace(); 173 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space"); 174 FESpace2Ceed(mesh_fes, ir, ceed, &mesh_basis, &mesh_restr); 175 176 CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &node_coords); 177 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, 178 mesh->GetNodes()->GetData()); 179 180 build_ctx.dim = mesh->Dimension(); 181 build_ctx.space_dim = mesh->SpaceDimension(); 182 183 CeedQFunctionCreateInterior(ceed, 1, 1, sizeof(CeedScalar), 184 (CeedEvalMode)(CEED_EVAL_GRAD|CEED_EVAL_WEIGHT), 185 CEED_EVAL_NONE, f_build_mass, 186 __FILE__":f_build_mass", &build_qfunc); 187 CeedQFunctionSetContext(build_qfunc, &build_ctx, sizeof(build_ctx)); 188 CeedOperatorCreate(ceed, mesh_restr, mesh_basis, build_qfunc, NULL, NULL, 189 &build_oper); 190 CeedOperatorGetQData(build_oper, &qdata); 191 CeedOperatorApply(build_oper, qdata, node_coords, NULL, 192 CEED_REQUEST_IMMEDIATE); 193 194 CeedQFunctionCreateInterior(ceed, 1, 1, sizeof(CeedScalar), 195 CEED_EVAL_INTERP, CEED_EVAL_INTERP, f_apply_mass, 196 __FILE__":f_apply_mass", &apply_qfunc); 197 CeedOperatorCreate(ceed, restr, basis, apply_qfunc, NULL, NULL, &oper); 198 199 CeedVectorCreate(ceed, fes->GetNDofs(), &u); 200 CeedVectorCreate(ceed, fes->GetNDofs(), &v); 201 } 202 203 /// Destructor 204 ~CeedMassOperator() { 205 CeedVectorDestroy(&v); 206 CeedVectorDestroy(&u); 207 CeedOperatorDestroy(&oper); 208 CeedQFunctionDestroy(&apply_qfunc); 209 // CeedVectorDestroy(&qdata); // qdata is owned by build_oper 210 CeedOperatorDestroy(&build_oper); 211 CeedQFunctionDestroy(&build_qfunc); 212 CeedVectorDestroy(&node_coords); 213 CeedElemRestrictionDestroy(&mesh_restr); 214 CeedBasisDestroy(&mesh_basis); 215 CeedElemRestrictionDestroy(&restr); 216 CeedBasisDestroy(&basis); 217 } 218 219 /// Operator action 220 virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const { 221 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_USE_POINTER, x.GetData()); 222 CeedVectorSetArray(v, CEED_MEM_HOST, CEED_USE_POINTER, y.GetData()); 223 CeedOperatorApply(oper, qdata, u, v, CEED_REQUEST_IMMEDIATE); 224 } 225 }; 226