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