xref: /libCEED/julia/LibCEED.jl/examples/common.jl (revision ffa5d67cac94379470c78ef400e8bd2c0655d3e7) 
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        ) :
90        nothing
91    restr_i =
92        form_strided ?
93        create_elem_restriction_strided(
94            c,
95            num_elem,
96            elem_qpts,
97            ncomp,
98            ncomp*elem_qpts*num_elem,
99            STRIDES_BACKEND,
100        ) :
101        nothing
102
103    return size, restr, restr_i
104end
105
106function set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
107    p = mesh_order
108    nd = p*nxyz .+ 1
109    num_elem::CeedInt = prod(nxyz)
110    scalar_size::CeedInt = prod(nd)
111
112    # The H1 basis uses Lobatto quadrature points as nodes.
113    nodes::Vector{CeedScalar} = lobatto_quadrature(p + 1) # nodes are in [-1,1]
114    nodes = 0.5 .+ 0.5*nodes
115
116    @witharray coords = mesh_coords begin
117        @inbounds @simd for gsnodes = 0:scalar_size-1
118            rnodes = gsnodes
119            for d = 1:dim
120                ndd = nd[d]
121                q = rnodes÷ndd
122                r = rnodes - ndd*q
123
124                d1d = r
125                coords[gsnodes+scalar_size*(d-1)+1] = ((d1d÷p) + nodes[d1d%p+1])/nxyz[d]
126                rnodes = q
127            end
128        end
129    end
130end
131