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 /// Mass operator example using MFEM 19 #include <ceed.h> 20 #include <mfem.hpp> 21 22 /// A structure used to pass additional data to f_build_mass 23 struct BuildContext { CeedInt dim, space_dim; }; 24 25 /// libCEED Q-function for building quadrature data for a mass operator 26 static int f_build_mass(void *ctx, CeedInt Q, 27 const CeedScalar *const *in, CeedScalar *const *out) { 28 // in[0] is Jacobians with shape [dim, nc=dim, Q] 29 // in[1] is quadrature weights, size (Q) 30 BuildContext *bc = (BuildContext *)ctx; 31 const CeedScalar *J = in[0], *qw = in[1]; 32 CeedScalar *rho = out[0]; 33 switch (bc->dim + 10*bc->space_dim) { 34 case 11: 35 for (CeedInt i=0; i<Q; i++) { 36 rho[i] = J[i] * qw[i]; 37 } 38 break; 39 case 22: 40 for (CeedInt i=0; i<Q; i++) { 41 // 0 2 42 // 1 3 43 rho[i] = (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * qw[i]; 44 } 45 break; 46 case 33: 47 for (CeedInt i=0; i<Q; i++) { 48 // 0 3 6 49 // 1 4 7 50 // 2 5 8 51 rho[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) - 52 J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) + 53 J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * qw[i]; 54 } 55 break; 56 default: 57 return CeedError(NULL, 1, "dim=%d, space_dim=%d is not supported", 58 bc->dim, bc->space_dim); 59 } 60 return 0; 61 } 62 63 /// libCEED Q-function for applying a mass operator 64 static int f_apply_mass(void *ctx, CeedInt Q, 65 const CeedScalar *const *in, CeedScalar *const *out) { 66 const CeedScalar *u = in[0], *w = in[1]; 67 CeedScalar *v = out[0]; 68 for (CeedInt i=0; i<Q; i++) { 69 v[i] = w[i] * u[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, restr_i, mesh_restr_i; 81 CeedQFunction apply_qfunc, build_qfunc; 82 CeedVector node_coords, rho; 83 CeedVector u, v; 84 85 BuildContext build_ctx; 86 87 static void FESpace2Ceed(const mfem::FiniteElementSpace *fes, 88 const mfem::IntegrationRule &ir, 89 Ceed ceed, CeedBasis *basis, 90 CeedElemRestriction *restr) { 91 mfem::Mesh *mesh = fes->GetMesh(); 92 const mfem::FiniteElement *fe = fes->GetFE(0); 93 const int order = fes->GetOrder(0); 94 mfem::Array<int> dof_map; 95 switch (mesh->Dimension()) { 96 case 1: { 97 const mfem::H1_SegmentElement *h1_fe = 98 dynamic_cast<const mfem::H1_SegmentElement *>(fe); 99 MFEM_VERIFY(h1_fe, "invalid FE"); 100 h1_fe->GetDofMap().Copy(dof_map); 101 break; 102 } 103 case 2: { 104 const mfem::H1_QuadrilateralElement *h1_fe = 105 dynamic_cast<const mfem::H1_QuadrilateralElement *>(fe); 106 MFEM_VERIFY(h1_fe, "invalid FE"); 107 h1_fe->GetDofMap().Copy(dof_map); 108 break; 109 } 110 case 3: { 111 const mfem::H1_HexahedronElement *h1_fe = 112 dynamic_cast<const mfem::H1_HexahedronElement *>(fe); 113 MFEM_VERIFY(h1_fe, "invalid FE"); 114 h1_fe->GetDofMap().Copy(dof_map); 115 break; 116 } 117 } 118 const mfem::FiniteElement *fe1d = 119 fes->FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT); 120 mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints()); 121 mfem::DenseMatrix grad1d(fe1d->GetDof(), ir.GetNPoints()); 122 mfem::Vector qref1d(ir.GetNPoints()), qweight1d(ir.GetNPoints()); 123 mfem::Vector shape_i(shape1d.Height()); 124 mfem::DenseMatrix grad_i(grad1d.Height(), 1); 125 const mfem::H1_SegmentElement *h1_fe1d = 126 dynamic_cast<const mfem::H1_SegmentElement *>(fe1d); 127 MFEM_VERIFY(h1_fe1d, "invalid FE"); 128 const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap(); 129 for (int i = 0; i < ir.GetNPoints(); i++) { 130 const mfem::IntegrationPoint &ip = ir.IntPoint(i); 131 qref1d(i) = ip.x; 132 qweight1d(i) = ip.weight; 133 fe1d->CalcShape(ip, shape_i); 134 fe1d->CalcDShape(ip, grad_i); 135 for (int j = 0; j < shape1d.Height(); j++) { 136 shape1d(j,i) = shape_i(dof_map_1d[j]); 137 grad1d(j,i) = grad_i(dof_map_1d[j],0); 138 } 139 } 140 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes->GetVDim(), order+1, 141 ir.GetNPoints(), shape1d.GetData(), 142 grad1d.GetData(), qref1d.GetData(), 143 qweight1d.GetData(), basis); 144 145 const mfem::Table &el_dof = fes->GetElementToDofTable(); 146 mfem::Array<int> tp_el_dof(el_dof.Size_of_connections()); 147 for (int i = 0; i < mesh->GetNE(); i++) { 148 const int el_offset = fe->GetDof()*i; 149 for (int j = 0; j < fe->GetDof(); j++) { 150 tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset]; 151 } 152 } 153 CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(), 154 fes->GetNDofs(), fes->GetVDim(), CEED_MEM_HOST, CEED_COPY_VALUES, 155 tp_el_dof.GetData(), restr); 156 } 157 158 public: 159 /// Constructor. Assumes @a fes is a scalar FE space. 160 CeedMassOperator(Ceed ceed, const mfem::FiniteElementSpace *fes) 161 : Operator(fes->GetNDofs()), 162 fes(fes) { 163 mfem::Mesh *mesh = fes->GetMesh(); 164 const int order = fes->GetOrder(0); 165 const int ir_order = 2*(order + 2) - 1; // <----- 166 const mfem::IntegrationRule &ir = 167 mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order); 168 CeedInt nqpts, nelem = mesh->GetNE(); 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 CeedBasisGetNumQuadraturePoints(basis, &nqpts); 176 177 CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts, 178 nqpts*nelem, 1, &restr_i); 179 CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts, 180 nqpts*nelem, 1, &mesh_restr_i); 181 182 CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &node_coords); 183 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, 184 mesh->GetNodes()->GetData()); 185 186 CeedVectorCreate(ceed, nelem*nqpts, &rho); 187 188 // Context data to be passed to the 'f_build_mass' Q-function. 189 build_ctx.dim = mesh->Dimension(); 190 build_ctx.space_dim = mesh->SpaceDimension(); 191 192 // Create the Q-function that builds the mass operator (i.e. computes its 193 // quadrature data) and set its context data. 194 CeedQFunctionCreateInterior(ceed, 1, f_build_mass, 195 __FILE__":f_build_mass", &build_qfunc); 196 CeedQFunctionAddInput(build_qfunc, "dx", mesh->SpaceDimension(), 197 CEED_EVAL_GRAD); 198 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT); 199 CeedQFunctionAddOutput(build_qfunc, "rho", 1, CEED_EVAL_NONE); 200 CeedQFunctionSetContext(build_qfunc, &build_ctx, sizeof(build_ctx)); 201 202 // Create the operator that builds the quadrature data for the mass operator. 203 CeedOperatorCreate(ceed, build_qfunc, NULL, NULL, &build_oper); 204 CeedOperatorSetField(build_oper, "dx", mesh_restr, CEED_NOTRANSPOSE, 205 mesh_basis, CEED_VECTOR_ACTIVE); 206 CeedOperatorSetField(build_oper, "weights", mesh_restr_i, CEED_NOTRANSPOSE, 207 mesh_basis, CEED_VECTOR_NONE); 208 CeedOperatorSetField(build_oper, "rho", restr_i, CEED_NOTRANSPOSE, 209 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 210 211 // Compute the quadrature data for the mass operator. 212 CeedOperatorApply(build_oper, node_coords, rho, 213 CEED_REQUEST_IMMEDIATE); 214 215 // Create the Q-function that defines the action of the mass operator. 216 CeedQFunctionCreateInterior(ceed, 1, f_apply_mass, 217 __FILE__":f_apply_mass", &apply_qfunc); 218 CeedQFunctionAddInput(apply_qfunc, "u", 1, CEED_EVAL_INTERP); 219 CeedQFunctionAddInput(apply_qfunc, "rho", 1, CEED_EVAL_NONE); 220 CeedQFunctionAddOutput(apply_qfunc, "v", 1, CEED_EVAL_INTERP); 221 222 // Create the mass operator. 223 CeedOperatorCreate(ceed, apply_qfunc, NULL, NULL, &oper); 224 CeedOperatorSetField(oper, "u", restr, CEED_NOTRANSPOSE, 225 basis, CEED_VECTOR_ACTIVE); 226 CeedOperatorSetField(oper, "rho", restr_i, CEED_NOTRANSPOSE, 227 CEED_BASIS_COLLOCATED, rho); 228 CeedOperatorSetField(oper, "v", restr, CEED_NOTRANSPOSE, 229 basis, CEED_VECTOR_ACTIVE); 230 231 CeedVectorCreate(ceed, fes->GetNDofs(), &u); 232 CeedVectorCreate(ceed, fes->GetNDofs(), &v); 233 } 234 235 /// Destructor 236 ~CeedMassOperator() { 237 CeedVectorDestroy(&u); 238 CeedVectorDestroy(&v); 239 CeedVectorDestroy(&rho); 240 CeedVectorDestroy(&node_coords); 241 CeedOperatorDestroy(&build_oper); 242 CeedQFunctionDestroy(&build_qfunc); 243 CeedOperatorDestroy(&oper); 244 CeedQFunctionDestroy(&apply_qfunc); 245 CeedBasisDestroy(&basis); 246 CeedBasisDestroy(&mesh_basis); 247 CeedElemRestrictionDestroy(&restr); 248 CeedElemRestrictionDestroy(&mesh_restr); 249 CeedElemRestrictionDestroy(&restr_i); 250 CeedElemRestrictionDestroy(&mesh_restr_i); 251 } 252 253 /// Operator action 254 virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const { 255 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_USE_POINTER, x.GetData()); 256 CeedVectorSetArray(v, CEED_MEM_HOST, CEED_USE_POINTER, y.GetData()); 257 CeedOperatorApply(oper, u, v, CEED_REQUEST_IMMEDIATE); 258 CeedVectorSyncArray(v, CEED_MEM_HOST); 259 } 260 }; 261