xref: /libCEED/julia/LibCEED.jl/examples/common.jl (revision daadeac6547c0bce0e170b8a41c931051f52e9a3) 
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{:restr}()
25const StridedOnly = FormRestrictionMode{:restr_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_restr = (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_restr
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    restr =
80        form_restr ?
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    restr_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, restr, restr_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