1function get_cartesian_mesh_size(dim, order, prob_size) 2 dims = zeros(CeedInt, dim) 3 # Use the approximate formula: 4 # prob_size ~ num_elem * order^dim 5 num_elem = div(prob_size, order^dim) 6 s = 0 # find s: num_elem/2 < 2^s <= num_elem 7 while num_elem > 1 8 num_elem = div(num_elem, 2) 9 s += 1 10 end 11 r = s%dim 12 for d = 1:dim 13 sd = div(s, dim) 14 if r > 0 15 sd += 1 16 r -= 1 17 end 18 dims[d] = 2^sd 19 end 20 dims 21end 22 23struct FormRestrictionMode{T} end 24const RestrictionOnly = FormRestrictionMode{:rstr}() 25const StridedOnly = FormRestrictionMode{:rstr_i}() 26const RestrictionAndStrided = FormRestrictionMode{:both}() 27 28function build_cartesian_restriction( 29 c::Ceed, 30 dim, 31 nxyz, 32 order, 33 ncomp, 34 num_qpts; 35 mode::Mode=RestrictionOnly, 36) where {Mode} 37 p::CeedInt = order 38 pp1::CeedInt = p + 1 39 nnodes::CeedInt = pp1^dim # number of scal. nodes per element 40 elem_qpts::CeedInt = num_qpts^dim # number of qpts per element 41 42 nd = CeedInt.(p*nxyz .+ 1) 43 num_elem::CeedInt = prod(nxyz) 44 scalar_size::CeedInt = prod(nd) 45 size::CeedInt = scalar_size*ncomp 46 47 form_rstr = (Mode() != StridedOnly) 48 form_strided = (Mode() != RestrictionOnly) 49 50 # elem: 0 1 n-1 51 # |---*-...-*---|---*-...-*---|- ... -|--...--| 52 # nnodes: 0 1 p-1 p p+1 2*p n*p 53 54 if form_rstr 55 el_nodes = zeros(CeedInt, num_elem*nnodes) 56 exyz = zeros(CeedInt, dim) 57 @inbounds @simd for e = 0:(num_elem-1) 58 re::CeedInt = e 59 for d::CeedInt = 1:dim 60 exyz[d] = re%nxyz[d] 61 re ÷= nxyz[d] 62 end 63 for lnodes::CeedInt = 0:(nnodes-1) 64 gnodes::CeedInt = 0 65 gnodes_stride::CeedInt = 1 66 rnodes::CeedInt = lnodes 67 for d::CeedInt = 1:dim 68 q = rnodes÷pp1 69 r = rnodes - pp1*q 70 gnodes::CeedInt += (exyz[d]*p + r)*gnodes_stride 71 gnodes_stride::CeedInt *= nd[d] 72 rnodes = q 73 end 74 el_nodes[e*nnodes+lnodes+1] = gnodes 75 end 76 end 77 end 78 79 rstr = 80 form_rstr ? 81 create_elem_restriction( 82 c, 83 num_elem, 84 nnodes, 85 ncomp, 86 scalar_size, 87 ncomp*scalar_size, 88 el_nodes, 89 ) : nothing 90 rstr_i = 91 form_strided ? 92 create_elem_restriction_strided( 93 c, 94 num_elem, 95 elem_qpts, 96 ncomp, 97 ncomp*elem_qpts*num_elem, 98 STRIDES_BACKEND, 99 ) : nothing 100 101 return size, rstr, rstr_i 102end 103 104function set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords) 105 p = mesh_order 106 nd = p*nxyz .+ 1 107 num_elem::CeedInt = prod(nxyz) 108 scalar_size::CeedInt = prod(nd) 109 110 # The H1 basis uses Lobatto quadrature points as nodes. 111 nodes::Vector{CeedScalar} = lobatto_quadrature(p + 1) # nodes are in [-1,1] 112 nodes = 0.5 .+ 0.5*nodes 113 114 @witharray coords = mesh_coords begin 115 @inbounds @simd for gsnodes = 0:scalar_size-1 116 rnodes = gsnodes 117 for d = 1:dim 118 ndd = nd[d] 119 q = rnodes÷ndd 120 r = rnodes - ndd*q 121 122 d1d = r 123 coords[gsnodes+scalar_size*(d-1)+1] = ((d1d÷p) + nodes[d1d%p+1])/nxyz[d] 124 rnodes = q 125 end 126 end 127 end 128end 129