xref: /libCEED/python/tests/test-3-basis.py (revision edf0491998c1d524f2f70fdd683669b8904cb3b6) 
1# Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors
2# All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3#
4# SPDX-License-Identifier: BSD-2-Clause
5#
6# This file is part of CEED:  http://github.com/ceed
7
8# @file
9# Test Ceed Basis functionality
10
11import os
12import math
13import libceed
14import numpy as np
15import buildmats as bm
16import check
17
18TOL = libceed.EPSILON * 256
19
20# -------------------------------------------------------------------------------
21# Utilities
22# -------------------------------------------------------------------------------
23
24
25def eval(dim, x):
26    result, center = 1, 0.1
27    for d in range(dim):
28        result *= math.tanh(x[d] - center)
29        center += 0.1
30    return result
31
32
33def feval(x1, x2):
34    return x1 * x1 + x2 * x2 + x1 * x2 + 1
35
36
37def dfeval(x1, x2):
38    return 2 * x1 + x2
39
40# -------------------------------------------------------------------------------
41# Test creation and distruction of a H1Lagrange basis
42# -------------------------------------------------------------------------------
43
44
45def test_300(ceed_resource, capsys):
46    ceed = libceed.Ceed(ceed_resource)
47
48    b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS_LOBATTO)
49    print(b)
50    del b
51
52    b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS)
53    print(b)
54    del b
55
56    # Only run this test in double precision
57    if libceed.lib.CEED_SCALAR_TYPE == libceed.SCALAR_FP64:
58        stdout, stderr, ref_stdout = check.output(capsys)
59        assert not stderr
60        assert stdout == ref_stdout
61
62# -------------------------------------------------------------------------------
63# Test GetNumNodes and GetNumQuadraturePoints for basis
64# -------------------------------------------------------------------------------
65
66
67def test_306(ceed_resource):
68    ceed = libceed.Ceed(ceed_resource)
69
70    b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO)
71
72    p = b.get_num_nodes()
73    q = b.get_num_quadrature_points()
74
75    assert p == 64
76    assert q == 125
77
78# -------------------------------------------------------------------------------
79# Test interpolation in multiple dimensions
80# -------------------------------------------------------------------------------
81
82
83def test_313(ceed_resource):
84    ceed = libceed.Ceed(ceed_resource)
85
86    for dim in range(1, 4):
87        Q = 10
88        Qdim = Q**dim
89        Xdim = 2**dim
90        x = np.empty(Xdim * dim, dtype=ceed.scalar_type())
91        uq = np.empty(Qdim, dtype=ceed.scalar_type())
92
93        for d in range(dim):
94            for i in range(Xdim):
95                x[d * Xdim + i] = 1 if (i %
96                                        (2**(dim - d))) // (2**(dim - d - 1)) else -1
97
98        X = ceed.Vector(Xdim * dim)
99        X.set_array(x, cmode=libceed.USE_POINTER)
100        Xq = ceed.Vector(Qdim * dim)
101        Xq.set_value(0)
102        U = ceed.Vector(Qdim)
103        U.set_value(0)
104        Uq = ceed.Vector(Qdim)
105
106        bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)
107        bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)
108
109        bxl.apply(1, libceed.EVAL_INTERP, X, Xq)
110
111        with Xq.array_read() as xq:
112            for i in range(Qdim):
113                xx = np.empty(dim, dtype=ceed.scalar_type())
114                for d in range(dim):
115                    xx[d] = xq[d * Qdim + i]
116                uq[i] = eval(dim, xx)
117
118        Uq.set_array(uq, cmode=libceed.USE_POINTER)
119
120        # This operation is the identity because the quadrature is collocated
121        bul.T.apply(1, libceed.EVAL_INTERP, Uq, U)
122
123        bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)
124        bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)
125
126        bxg.apply(1, libceed.EVAL_INTERP, X, Xq)
127        bug.apply(1, libceed.EVAL_INTERP, U, Uq)
128
129        with Xq.array_read() as xq, Uq.array_read() as u:
130            for i in range(Qdim):
131                xx = np.empty(dim, dtype=ceed.scalar_type())
132                for d in range(dim):
133                    xx[d] = xq[d * Qdim + i]
134                fx = eval(dim, xx)
135                assert math.fabs(u[i] - fx) < 1E-4
136
137# -------------------------------------------------------------------------------
138# Test grad in multiple dimensions
139# -------------------------------------------------------------------------------
140
141
142def test_314(ceed_resource):
143    ceed = libceed.Ceed(ceed_resource)
144
145    for dim in range(1, 4):
146        P, Q = 8, 10
147        Pdim = P**dim
148        Qdim = Q**dim
149        Xdim = 2**dim
150        sum1 = sum2 = 0
151        x = np.empty(Xdim * dim, dtype=ceed.scalar_type())
152        u = np.empty(Pdim, dtype=ceed.scalar_type())
153
154        for d in range(dim):
155            for i in range(Xdim):
156                x[d * Xdim + i] = 1 if (i %
157                                        (2**(dim - d))) // (2**(dim - d - 1)) else -1
158
159        X = ceed.Vector(Xdim * dim)
160        X.set_array(x, cmode=libceed.USE_POINTER)
161        Xq = ceed.Vector(Pdim * dim)
162        Xq.set_value(0)
163        U = ceed.Vector(Pdim)
164        Uq = ceed.Vector(Qdim * dim)
165        Uq.set_value(0)
166        Ones = ceed.Vector(Qdim * dim)
167        Ones.set_value(1)
168        Gtposeones = ceed.Vector(Pdim)
169        Gtposeones.set_value(0)
170
171        # Get function values at quadrature points
172        bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)
173        bxl.apply(1, libceed.EVAL_INTERP, X, Xq)
174
175        with Xq.array_read() as xq:
176            for i in range(Pdim):
177                xx = np.empty(dim, dtype=ceed.scalar_type())
178                for d in range(dim):
179                    xx[d] = xq[d * Pdim + i]
180                u[i] = eval(dim, xx)
181
182        U.set_array(u, cmode=libceed.USE_POINTER)
183
184        # Calculate G u at quadrature points, G' * 1 at dofs
185        bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)
186        bug.apply(1, libceed.EVAL_GRAD, U, Uq)
187        bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones)
188
189        # Check if 1' * G * u = u' * (G' * 1)
190        with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq:
191            for i in range(Pdim):
192                sum1 += gtposeones[i] * u[i]
193            for i in range(dim * Qdim):
194                sum2 += uq[i]
195
196        assert math.fabs(sum1 - sum2) < 10. * TOL
197
198# -------------------------------------------------------------------------------
199# Test creation and destruction of a 2D Simplex non-tensor H1 basis
200# -------------------------------------------------------------------------------
201
202
203def test_320(ceed_resource):
204    ceed = libceed.Ceed(ceed_resource)
205
206    P, Q, dim = 6, 4, 2
207
208    in_array = np.empty(P, dtype=ceed.scalar_type())
209    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
210    qweight = np.empty(Q, dtype=ceed.scalar_type())
211
212    interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[
213        libceed.lib.CEED_SCALAR_TYPE])
214
215    b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight)
216
217    print(b)
218    del b
219
220# -------------------------------------------------------------------------------
221# Test integration with a 2D Simplex non-tensor H1 basis
222# -------------------------------------------------------------------------------
223
224
225def test_322(ceed_resource):
226    ceed = libceed.Ceed(ceed_resource)
227
228    P, Q, dim = 6, 4, 2
229
230    xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0.,
231                   0., 0.5, 0.5, 1.], dtype=ceed.scalar_type())
232    in_array = np.empty(P, dtype=ceed.scalar_type())
233    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
234    qweight = np.empty(Q, dtype=ceed.scalar_type())
235
236    interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[
237        libceed.lib.CEED_SCALAR_TYPE])
238
239    b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight)
240
241    # Interpolate function to quadrature points
242    for i in range(P):
243        in_array[i] = feval(xr[0 * P + i], xr[1 * P + i])
244
245    in_vec = ceed.Vector(P)
246    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
247    out_vec = ceed.Vector(Q)
248    out_vec.set_value(0)
249    weights_vec = ceed.Vector(Q)
250    weights_vec.set_value(0)
251
252    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
253    b.apply(1, libceed.EVAL_WEIGHT, libceed.VECTOR_NONE, weights_vec)
254
255    # Check values at quadrature points
256    with out_vec.array_read() as out_array, weights_vec.array_read() as weights_array:
257        sum = 0
258        for i in range(Q):
259            sum += out_array[i] * weights_array[i]
260        assert math.fabs(sum - 17. / 24.) < 10. * TOL
261
262# -------------------------------------------------------------------------------
263# Test grad with a 2D Simplex non-tensor H1 basis
264# -------------------------------------------------------------------------------
265
266
267def test_323(ceed_resource):
268    ceed = libceed.Ceed(ceed_resource)
269
270    P, Q, dim = 6, 4, 2
271
272    xq = np.array([0.2, 0.6, 1. / 3., 0.2, 0.2, 0.2,
273                   1. / 3., 0.6], dtype=ceed.scalar_type())
274    xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0.,
275                   0., 0.5, 0.5, 1.], dtype=ceed.scalar_type())
276    in_array = np.empty(P, dtype=ceed.scalar_type())
277    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
278    qweight = np.empty(Q, dtype=ceed.scalar_type())
279
280    interp, grad = bm.buildmats(qref, qweight, libceed.scalar_types[
281        libceed.lib.CEED_SCALAR_TYPE])
282
283    b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight)
284
285    # Interpolate function to quadrature points
286    for i in range(P):
287        in_array[i] = feval(xr[0 * P + i], xr[1 * P + i])
288
289    in_vec = ceed.Vector(P)
290    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
291    out_vec = ceed.Vector(Q * dim)
292    out_vec.set_value(0)
293
294    b.apply(1, libceed.EVAL_GRAD, in_vec, out_vec)
295
296    # Check values at quadrature points
297    with out_vec.array_read() as out_array:
298        for i in range(Q):
299            value = dfeval(xq[0 * Q + i], xq[1 * Q + i])
300            assert math.fabs(out_array[0 * Q + i] - value) < 10. * TOL
301
302            value = dfeval(xq[1 * Q + i], xq[0 * Q + i])
303            assert math.fabs(out_array[1 * Q + i] - value) < 10. * TOL
304
305# -------------------------------------------------------------------------------
306# Test interpolation with a 2D Quad non-tensor H(div) basis
307# -------------------------------------------------------------------------------
308
309
310def test_330(ceed_resource):
311    ceed = libceed.Ceed(ceed_resource)
312
313    P, Q, dim = 4, 4, 2
314
315    in_array = np.ones(P, dtype=ceed.scalar_type())
316    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
317    qweight = np.empty(Q, dtype=ceed.scalar_type())
318
319    interp, div = bm.buildmatshdiv(qref, qweight, libceed.scalar_types[
320        libceed.lib.CEED_SCALAR_TYPE])
321
322    b = ceed.BasisHdiv(libceed.QUAD, 1, P, Q, interp, div, qref, qweight)
323
324    # Interpolate function to quadrature points
325    in_vec = ceed.Vector(P)
326    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
327    out_vec = ceed.Vector(Q * dim)
328    out_vec.set_value(0)
329
330    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
331
332    # Check values at quadrature points
333    with out_vec.array_read() as out_array:
334        for i in range(Q):
335            assert math.fabs(out_array[0 * Q + i] + 1.) < 10. * TOL
336            assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL
337
338# -------------------------------------------------------------------------------
339# Test interpolation with a 2D Simplex non-tensor H(curl) basis
340# -------------------------------------------------------------------------------
341
342
343def test_340(ceed_resource):
344    ceed = libceed.Ceed(ceed_resource)
345
346    P, Q, dim = 3, 4, 2
347
348    in_array = np.empty(P, dtype=ceed.scalar_type())
349    qref = np.empty(dim * Q, dtype=ceed.scalar_type())
350    qweight = np.empty(Q, dtype=ceed.scalar_type())
351
352    interp, curl = bm.buildmatshcurl(qref, qweight, libceed.scalar_types[
353        libceed.lib.CEED_SCALAR_TYPE])
354
355    b = ceed.BasisHcurl(libceed.TRIANGLE, 1, P, Q, interp, curl, qref, qweight)
356
357    # Interpolate function to quadrature points
358    in_array[0] = 1.
359    in_array[1] = 2.
360    in_array[2] = 1.
361
362    in_vec = ceed.Vector(P)
363    in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
364    out_vec = ceed.Vector(Q * dim)
365    out_vec.set_value(0)
366
367    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
368
369    # Check values at quadrature points
370    with out_vec.array_read() as out_array:
371        for i in range(Q):
372            assert math.fabs(out_array[0 * Q + i] - 1.) < 10. * TOL
373
374    # Interpolate function to quadrature points
375    in_array[0] = -1.
376    in_array[1] = 1.
377    in_array[2] = 2.
378
379    out_vec.set_value(0)
380
381    b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
382
383    # Check values at quadrature points
384    with out_vec.array_read() as out_array:
385        for i in range(Q):
386            assert math.fabs(out_array[1 * Q + i] - 1.) < 10. * TOL
387
388# -------------------------------------------------------------------------------
389