1using LibCEED.C, Printf, UnsafeArrays 2using LibCEED.C: CeedInt, CeedScalar 3 4# A structure used to pass additional data to f_build_mass 5mutable struct BuildContextC 6 dim::CeedInt 7 space_dim::CeedInt 8end 9 10# libCEED Q-function for building quadrature data for a mass operator 11function f_build_mass_c( 12 ctx_ptr::Ptr{Cvoid}, 13 Q::CeedInt, 14 in_ptr::Ptr{Ptr{CeedScalar}}, 15 out_ptr::Ptr{Ptr{CeedScalar}}, 16) 17 # in[0] is Jacobians with shape [dim, nc=dim, Q] 18 # in[1] is quadrature weights, size (Q) 19 ctx = unsafe_load(Ptr{BuildContextC}(ctx_ptr)) 20 J = UnsafeArray(unsafe_load(in_ptr, 1), (Int(Q), Int(ctx.dim^2))) 21 w = UnsafeArray(unsafe_load(in_ptr, 2), (Int(Q),)) 22 qdata = UnsafeArray(unsafe_load(out_ptr, 1), (Int(Q),)) 23 if ctx.dim == 1 24 @inbounds @simd for i = 1:Q 25 qdata[i] = J[i]*w[i] 26 end 27 elseif ctx.dim == 2 28 @inbounds @simd for i = 1:Q 29 qdata[i] = (J[i, 1]*J[i, 4] - J[i, 2]*J[i, 3])*w[i] 30 end 31 elseif ctx.dim == 3 32 @inbounds @simd for i = 1:Q 33 qdata[i] = 34 ( 35 J[i, 1]*(J[i, 5]*J[i, 9] - J[i, 6]*J[i, 8]) - 36 J[i, 2]*(J[i, 4]*J[i, 9] - J[i, 6]*J[i, 7]) + 37 J[i, 3]*(J[i, 4]*J[i, 8] - J[i, 5]*J[i, 7]) 38 )*w[i] 39 end 40 else 41 error("Bad dimension") 42 end 43 return CeedInt(0) 44end 45 46# libCEED Q-function for applying a mass operator 47function f_apply_mass_c( 48 ctx, 49 Q::CeedInt, 50 in_ptr::Ptr{Ptr{CeedScalar}}, 51 out_ptr::Ptr{Ptr{CeedScalar}}, 52) 53 u = UnsafeArray(unsafe_load(in_ptr, 1), (Int(Q),)) 54 qdata = UnsafeArray(unsafe_load(in_ptr, 2), (Int(Q),)) 55 v = UnsafeArray(unsafe_load(out_ptr, 1), (Int(Q),)) 56 @inbounds @simd for i = 1:Q 57 v[i] = qdata[i]*u[i] 58 end 59 return CeedInt(0) 60end 61 62function get_cartesian_mesh_size_c(dim, order, prob_size) 63 dims = zeros(Int, dim) 64 # Use the approximate formula: 65 # prob_size ~ num_elem * order^dim 66 num_elem = div(prob_size, order^dim) 67 s = 0 # find s: num_elem/2 < 2^s <= num_elem 68 while num_elem > 1 69 num_elem = div(num_elem, 2) 70 s += 1 71 end 72 r = s%dim 73 for d = 1:dim 74 sd = div(s, dim) 75 if r > 0 76 sd += 1 77 r -= 1 78 end 79 dims[d] = 2^sd 80 end 81 dims 82end 83 84function build_cartesian_restriction_c( 85 ceed, 86 dim, 87 nxyz, 88 order, 89 ncomp, 90 num_qpts; 91 form_strided=false, 92) 93 p = order 94 pp1 = p + 1 95 nnodes = pp1^dim # number of scal. nodes per element 96 elem_qpts = num_qpts^dim # number of qpts per element 97 num_elem = 1 98 scalar_size = 1 99 100 nd = p*nxyz .+ 1 101 num_elem = prod(nxyz) 102 scalar_size = prod(nd) 103 size = scalar_size*ncomp 104 105 # elem: 0 1 n-1 106 # |---*-...-*---|---*-...-*---|- ... -|--...--| 107 # nnodes: 0 1 p-1 p p+1 2*p n*p 108 109 el_nodes = zeros(C.CeedInt, num_elem*nnodes) 110 exyz = zeros(Int, dim) 111 @inbounds for e = 0:(num_elem-1) 112 re = e 113 for d = 1:dim 114 exyz[d] = re%nxyz[d] 115 re = div(re, nxyz[d]) 116 end 117 for lnodes = 0:(nnodes-1) 118 gnodes = 0 119 gnodes_stride = 1 120 rnodes = lnodes 121 for d = 1:dim 122 gnodes += (exyz[d]*p + rnodes%pp1)*gnodes_stride 123 gnodes_stride *= nd[d] 124 rnodes = div(rnodes, pp1) 125 end 126 el_nodes[e*nnodes+lnodes+1] = gnodes 127 end 128 end 129 130 rstr = Ref{C.CeedElemRestriction}() 131 C.CeedElemRestrictionCreate( 132 ceed[], 133 num_elem, 134 nnodes, 135 ncomp, 136 scalar_size, 137 ncomp*scalar_size, 138 C.CEED_MEM_HOST, 139 C.CEED_COPY_VALUES, 140 el_nodes, 141 rstr, 142 ) 143 if form_strided 144 rstr_i = Ref{C.CeedElemRestriction}() 145 err = C.CeedElemRestrictionCreateStrided( 146 ceed[], 147 num_elem, 148 elem_qpts, 149 ncomp, 150 ncomp*elem_qpts*num_elem, 151 C.CEED_STRIDES_BACKEND[], 152 rstr_i, 153 ) 154 return size, rstr, rstr_i 155 else 156 return size, rstr 157 end 158end 159 160function set_cartesian_mesh_coords_c(dim, nxyz, mesh_order, mesh_coords) 161 p = mesh_order 162 nd = p*nxyz .+ 1 163 num_elem = prod(nxyz) 164 scalar_size = prod(nd) 165 166 coords_ref = Ref{Ptr{C.CeedScalar}}() 167 C.CeedVectorGetArray(mesh_coords[], C.CEED_MEM_HOST, coords_ref) 168 coords = unsafe_wrap(Array, coords_ref[], scalar_size*dim) 169 170 nodes = zeros(C.CeedScalar, p + 1) 171 # The H1 basis uses Lobatto quadrature points as nodes. 172 C.CeedLobattoQuadrature(p + 1, nodes, C_NULL) # nodes are in [-1,1] 173 nodes = 0.5 .+ 0.5*nodes 174 for gsnodes = 0:(scalar_size-1) 175 rnodes = gsnodes 176 for d = 1:dim 177 d1d = rnodes%nd[d] 178 coords[gsnodes+scalar_size*(d-1)+1] = (div(d1d, p) + nodes[d1d%p+1])/nxyz[d] 179 rnodes = div(rnodes, nd[d]) 180 end 181 end 182 C.CeedVectorRestoreArray(mesh_coords[], coords_ref) 183end 184 185function transform_mesh_coords_c(dim, mesh_size, mesh_coords) 186 coords_ref = Ref{Ptr{C.CeedScalar}}() 187 C.CeedVectorGetArray(mesh_coords[], C.CEED_MEM_HOST, coords_ref) 188 coords = unsafe_wrap(Array, coords_ref[], mesh_size) 189 190 if dim == 1 191 for i = 1:mesh_size 192 # map [0,1] to [0,1] varying the mesh density 193 coords[i] = 0.5 + 1.0/sqrt(3.0)*sin((2.0/3.0)*pi*(coords[i] - 0.5)) 194 end 195 exact_volume = 1 196 else 197 num_nodes = div(mesh_size, dim) 198 for i = 1:num_nodes 199 # map (x,y) from [0,1]x[0,1] to the quarter annulus with polar 200 # coordinates, (r,phi) in [1,2]x[0,pi/2] with area = 3/4*pi 201 u = coords[i] 202 v = coords[i+num_nodes] 203 u = 1.0 + u 204 v = pi/2*v 205 coords[i] = u*cos(v) 206 coords[i+num_nodes] = u*sin(v) 207 end 208 exact_volume = 3.0/4.0*pi 209 end 210 211 C.CeedVectorRestoreArray(mesh_coords[], coords_ref) 212 return exact_volume 213end 214 215function run_ex1_c(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size) 216 ncompx = dim 217 prob_size < 0 && (prob_size = 256*1024) 218 219 gallery = false 220 221 ceed = Ref{C.Ceed}() 222 C.CeedInit(ceed_spec, ceed) 223 224 mesh_basis = Ref{C.CeedBasis}() 225 sol_basis = Ref{C.CeedBasis}() 226 C.CeedBasisCreateTensorH1Lagrange( 227 ceed[], 228 dim, 229 ncompx, 230 mesh_order + 1, 231 num_qpts, 232 C.CEED_GAUSS, 233 mesh_basis, 234 ) 235 C.CeedBasisCreateTensorH1Lagrange( 236 ceed[], 237 dim, 238 1, 239 sol_order + 1, 240 num_qpts, 241 C.CEED_GAUSS, 242 sol_basis, 243 ) 244 245 # Determine the mesh size based on the given approximate problem size. 246 nxyz = get_cartesian_mesh_size_c(dim, sol_order, prob_size) 247 println("Mesh size: ", nxyz) 248 249 # Build CeedElemRestriction objects describing the mesh and solution discrete 250 # representations. 251 mesh_size, mesh_rstr = 252 build_cartesian_restriction_c(ceed, dim, nxyz, mesh_order, ncompx, num_qpts) 253 sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction_c( 254 ceed, 255 dim, 256 nxyz, 257 sol_order, 258 1, 259 num_qpts, 260 form_strided=true, 261 ) 262 println("Number of mesh nodes : ", div(mesh_size, dim)) 263 println("Number of solution nodes : ", sol_size) 264 265 # Create a C.CeedVector with the mesh coordinates. 266 mesh_coords = Ref{C.CeedVector}() 267 C.CeedVectorCreate(ceed[], mesh_size, mesh_coords) 268 set_cartesian_mesh_coords_c(dim, nxyz, mesh_order, mesh_coords) 269 # Apply a transformation to the mesh. 270 exact_vol = transform_mesh_coords_c(dim, mesh_size, mesh_coords) 271 272 # Create the Q-function that builds the mass operator (i.e. computes its 273 # quadrature data) and set its context data. 274 build_qfunc = Ref{C.CeedQFunction}() 275 276 build_ctx = BuildContextC(dim, dim) 277 qf_ctx = Ref{C.CeedQFunctionContext}() 278 C.CeedQFunctionContextCreate(ceed[], qf_ctx) 279 C.CeedQFunctionContextSetData( 280 qf_ctx[], 281 C.CEED_MEM_HOST, 282 C.CEED_USE_POINTER, 283 sizeof(build_ctx), 284 pointer_from_objref(build_ctx), 285 ) 286 287 if !gallery 288 qf_build_mass = @cfunction( 289 f_build_mass_c, 290 C.CeedInt, 291 (Ptr{Cvoid}, C.CeedInt, Ptr{Ptr{C.CeedScalar}}, Ptr{Ptr{C.CeedScalar}}) 292 ) 293 # This creates the QFunction directly. 294 C.CeedQFunctionCreateInterior(ceed[], 1, qf_build_mass, "julia", build_qfunc) 295 C.CeedQFunctionAddInput(build_qfunc[], "dx", ncompx*dim, C.CEED_EVAL_GRAD) 296 C.CeedQFunctionAddInput(build_qfunc[], "weights", 1, C.CEED_EVAL_WEIGHT) 297 C.CeedQFunctionAddOutput(build_qfunc[], "qdata", 1, C.CEED_EVAL_NONE) 298 C.CeedQFunctionSetContext(build_qfunc[], qf_ctx[]) 299 else 300 # This creates the QFunction via the gallery. 301 name = "Mass$(dim)DBuild" 302 C.CeedQFunctionCreateInteriorByName(ceed[], name, build_qfunc) 303 end 304 305 # Create the operator that builds the quadrature data for the mass operator. 306 build_oper = Ref{C.CeedOperator}() 307 C.CeedOperatorCreate( 308 ceed[], 309 build_qfunc[], 310 C.CEED_QFUNCTION_NONE[], 311 C.CEED_QFUNCTION_NONE[], 312 build_oper, 313 ) 314 C.CeedOperatorSetField( 315 build_oper[], 316 "dx", 317 mesh_rstr[], 318 mesh_basis[], 319 C.CEED_VECTOR_ACTIVE[], 320 ) 321 C.CeedOperatorSetField( 322 build_oper[], 323 "weights", 324 C.CEED_ELEMRESTRICTION_NONE[], 325 mesh_basis[], 326 C.CEED_VECTOR_NONE[], 327 ) 328 C.CeedOperatorSetField( 329 build_oper[], 330 "qdata", 331 sol_rstr_i[], 332 C.CEED_BASIS_NONE[], 333 C.CEED_VECTOR_ACTIVE[], 334 ) 335 336 # Compute the quadrature data for the mass operator. 337 qdata = Ref{C.CeedVector}() 338 elem_qpts = num_qpts^dim 339 num_elem = prod(nxyz) 340 C.CeedVectorCreate(ceed[], num_elem*elem_qpts, qdata) 341 342 print("Computing the quadrature data for the mass operator ...") 343 flush(stdout) 344 GC.@preserve build_ctx C.CeedOperatorApply( 345 build_oper[], 346 mesh_coords[], 347 qdata[], 348 C.CEED_REQUEST_IMMEDIATE[], 349 ) 350 println(" done.") 351 352 # Create the Q-function that defines the action of the mass operator. 353 apply_qfunc = Ref{C.CeedQFunction}() 354 if !gallery 355 qf_apply_mass = @cfunction( 356 f_apply_mass_c, 357 C.CeedInt, 358 (Ptr{Cvoid}, C.CeedInt, Ptr{Ptr{C.CeedScalar}}, Ptr{Ptr{C.CeedScalar}}) 359 ) 360 # This creates the QFunction directly. 361 C.CeedQFunctionCreateInterior(ceed[], 1, qf_apply_mass, "julia", apply_qfunc) 362 C.CeedQFunctionAddInput(apply_qfunc[], "u", 1, C.CEED_EVAL_INTERP) 363 C.CeedQFunctionAddInput(apply_qfunc[], "qdata", 1, C.CEED_EVAL_NONE) 364 C.CeedQFunctionAddOutput(apply_qfunc[], "v", 1, C.CEED_EVAL_INTERP) 365 else 366 # This creates the QFunction via the gallery. 367 C.CeedQFunctionCreateInteriorByName(ceed[], "MassApply", apply_qfunc) 368 end 369 370 # Create the mass operator. 371 oper = Ref{C.CeedOperator}() 372 C.CeedOperatorCreate( 373 ceed[], 374 apply_qfunc[], 375 C.CEED_QFUNCTION_NONE[], 376 C.CEED_QFUNCTION_NONE[], 377 oper, 378 ) 379 C.CeedOperatorSetField(oper[], "u", sol_rstr[], sol_basis[], C.CEED_VECTOR_ACTIVE[]) 380 C.CeedOperatorSetField(oper[], "qdata", sol_rstr_i[], C.CEED_BASIS_NONE[], qdata[]) 381 C.CeedOperatorSetField(oper[], "v", sol_rstr[], sol_basis[], C.CEED_VECTOR_ACTIVE[]) 382 383 # Compute the mesh volume using the mass operator: vol = 1^T \cdot M \cdot 1 384 print("Computing the mesh volume using the formula: vol = 1^T.M.1 ...") 385 flush(stdout) 386 # Create auxiliary solution-size vectors. 387 u = Ref{C.CeedVector}() 388 v = Ref{C.CeedVector}() 389 C.CeedVectorCreate(ceed[], sol_size, u) 390 C.CeedVectorCreate(ceed[], sol_size, v) 391 392 # Initialize 'u' with ones. 393 C.CeedVectorSetValue(u[], 1.0) 394 395 # Apply the mass operator: 'u' -> 'v'. 396 C.CeedOperatorApply(oper[], u[], v[], C.CEED_REQUEST_IMMEDIATE[]) 397 398 # Compute and print the sum of the entries of 'v' giving the mesh volume. 399 v_host_ref = Ref{Ptr{C.CeedScalar}}() 400 C.CeedVectorGetArrayRead(v[], C.CEED_MEM_HOST, v_host_ref) 401 v_host = unsafe_wrap(Array, v_host_ref[], sol_size) 402 vol = sum(v_host) 403 C.CeedVectorRestoreArrayRead(v[], v_host_ref) 404 405 println(" done.") 406 @printf("Exact mesh volume : % .14g\n", exact_vol) 407 @printf("Computed mesh volume : % .14g\n", vol) 408 @printf("Volume error : % .14g\n", vol - exact_vol) 409 410 # Free dynamically allocated memory. 411 C.CeedVectorDestroy(u) 412 C.CeedVectorDestroy(v) 413 C.CeedVectorDestroy(qdata) 414 C.CeedVectorDestroy(mesh_coords) 415 C.CeedOperatorDestroy(oper) 416 C.CeedQFunctionDestroy(apply_qfunc) 417 C.CeedOperatorDestroy(build_oper) 418 C.CeedQFunctionDestroy(build_qfunc) 419 C.CeedElemRestrictionDestroy(sol_rstr) 420 C.CeedElemRestrictionDestroy(mesh_rstr) 421 C.CeedElemRestrictionDestroy(sol_rstr_i) 422 C.CeedBasisDestroy(sol_basis) 423 C.CeedBasisDestroy(mesh_basis) 424 C.CeedDestroy(ceed) 425end 426 427run_ex1_c( 428 ceed_spec="/cpu/self", 429 dim=3, 430 mesh_order=4, 431 sol_order=4, 432 num_qpts=4 + 2, 433 prob_size=-1, 434) 435