xref: /libCEED/python/ceed_basis.py (revision a550dd668ad5fdaeaade4afb48509e4b41c51584) 
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
8from _ceed_cffi import ffi, lib
9import tempfile
10import numpy as np
11from abc import ABC
12from .ceed_constants import TRANSPOSE, NOTRANSPOSE, scalar_types
13
14# ------------------------------------------------------------------------------
15
16
17class Basis(ABC):
18    """Ceed Basis: finite element basis objects."""
19
20    # Representation
21    def __repr__(self):
22        return "<CeedBasis instance at " + hex(id(self)) + ">"
23
24    # String conversion for print() to stdout
25    def __str__(self):
26        """View a Basis via print()."""
27
28        # libCEED call
29        with tempfile.NamedTemporaryFile() as key_file:
30            with open(key_file.name, 'r+') as stream_file:
31                stream = ffi.cast("FILE *", stream_file)
32
33                err_code = lib.CeedBasisView(self._pointer[0], stream)
34                self._ceed._check_error(err_code)
35
36                stream_file.seek(0)
37                out_string = stream_file.read()
38
39        return out_string
40
41    # Apply Basis
42    def apply(self, nelem, emode, u, v, tmode=NOTRANSPOSE):
43        """Apply basis evaluation from nodes to quadrature points or vice versa.
44
45           Args:
46             nelem: the number of elements to apply the basis evaluation to;
47                      the backend will specify the ordering in a
48                      BlockedElemRestriction
49             emode: basis evaluation mode
50             u: input vector
51             v: output vector
52             **tmode: CEED_NOTRANSPOSE to evaluate from nodes to quadrature
53                        points, CEED_TRANSPOSE to apply the transpose, mapping
54                        from quadrature points to nodes; default CEED_NOTRANSPOSE"""
55
56        # libCEED call
57        err_code = lib.CeedBasisApply(self._pointer[0], nelem, tmode, emode,
58                                      u._pointer[0], v._pointer[0])
59        self._ceed._check_error(err_code)
60
61    # Transpose a Basis
62    @property
63    def T(self):
64        """Transpose a Basis."""
65
66        return TransposeBasis(self)
67
68    # Transpose a Basis
69    @property
70    def transpose(self):
71        """Transpose a Basis."""
72
73        return TransposeBasis(self)
74
75    # Get number of nodes
76    def get_num_nodes(self):
77        """Get total number of nodes (in dim dimensions) of a Basis.
78
79           Returns:
80             num_nodes: total number of nodes"""
81
82        # Setup argument
83        p_pointer = ffi.new("CeedInt *")
84
85        # libCEED call
86        err_code = lib.CeedBasisGetNumNodes(self._pointer[0], p_pointer)
87        self._ceed._check_error(err_code)
88
89        return p_pointer[0]
90
91    # Get number of quadrature points
92    def get_num_quadrature_points(self):
93        """Get total number of quadrature points (in dim dimensions) of a Basis.
94
95           Returns:
96             num_qpts: total number of quadrature points"""
97
98        # Setup argument
99        q_pointer = ffi.new("CeedInt *")
100
101        # libCEED call
102        err_code = lib.CeedBasisGetNumQuadraturePoints(
103            self._pointer[0], q_pointer)
104        self._ceed._check_error(err_code)
105
106        return q_pointer[0]
107
108    # Destructor
109    def __del__(self):
110        # libCEED call
111        err_code = lib.CeedBasisDestroy(self._pointer)
112        self._ceed._check_error(err_code)
113
114# ------------------------------------------------------------------------------
115
116
117class BasisTensorH1(Basis):
118    """Ceed Tensor H1 Basis: finite element tensor-product basis objects for
119         H^1 discretizations."""
120
121    # Constructor
122    def __init__(self, ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d,
123                 qref1d, qweight1d):
124
125        # Setup arguments
126        self._pointer = ffi.new("CeedBasis *")
127
128        self._ceed = ceed
129
130        interp1d_pointer = ffi.new("CeedScalar *")
131        interp1d_pointer = ffi.cast(
132            "CeedScalar *",
133            interp1d.__array_interface__['data'][0])
134
135        grad1d_pointer = ffi.new("CeedScalar *")
136        grad1d_pointer = ffi.cast(
137            "CeedScalar *",
138            grad1d.__array_interface__['data'][0])
139
140        qref1d_pointer = ffi.new("CeedScalar *")
141        qref1d_pointer = ffi.cast(
142            "CeedScalar *",
143            qref1d.__array_interface__['data'][0])
144
145        qweight1d_pointer = ffi.new("CeedScalar *")
146        qweight1d_pointer = ffi.cast(
147            "CeedScalar *",
148            qweight1d.__array_interface__['data'][0])
149
150        # libCEED call
151        err_code = lib.CeedBasisCreateTensorH1(self._ceed._pointer[0], dim, ncomp,
152                                               P1d, Q1d, interp1d_pointer,
153                                               grad1d_pointer, qref1d_pointer,
154                                               qweight1d_pointer, self._pointer)
155        self._ceed._check_error(err_code)
156
157    # Get 1D interpolation matrix
158    def get_interp_1d(self):
159        """Return 1D interpolation matrix of a tensor product Basis.
160
161           Returns:
162             *array: Numpy array"""
163
164        # Compute the length of the array
165        nnodes_pointer = ffi.new("CeedInt *")
166        lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer)
167        nqpts_pointer = ffi.new("CeedInt *")
168        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
169        length = nnodes_pointer[0] * nqpts_pointer[0]
170
171        # Setup the pointer's pointer
172        array_pointer = ffi.new("CeedScalar **")
173
174        # libCEED call
175        lib.CeedBasisGetInterp1D(self._pointer[0], array_pointer)
176
177        # Return array created from buffer
178        # Create buffer object from returned pointer
179        buff = ffi.buffer(
180            array_pointer[0],
181            ffi.sizeof("CeedScalar") *
182            length)
183        # return read only Numpy array
184        ret = np.frombuffer(buff, dtype=scalar_types[lib.CEED_SCALAR_TYPE])
185        ret.flags['WRITEABLE'] = False
186        return ret
187
188    # Get 1D gradient matrix
189    def get_grad_1d(self):
190        """Return 1D gradient matrix of a tensor product Basis.
191
192           Returns:
193             *array: Numpy array"""
194
195        # Compute the length of the array
196        nnodes_pointer = ffi.new("CeedInt *")
197        lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer)
198        nqpts_pointer = ffi.new("CeedInt *")
199        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
200        length = nnodes_pointer[0] * nqpts_pointer[0]
201
202        # Setup the pointer's pointer
203        array_pointer = ffi.new("CeedScalar **")
204
205        # libCEED call
206        lib.CeedBasisGetGrad1D(self._pointer[0], array_pointer)
207
208        # Return array created from buffer
209        # Create buffer object from returned pointer
210        buff = ffi.buffer(
211            array_pointer[0],
212            ffi.sizeof("CeedScalar") *
213            length)
214        # return read only Numpy array
215        ret = np.frombuffer(buff, dtype="float64")
216        ret.flags['WRITEABLE'] = False
217        return ret
218
219    # Get 1D quadrature weights matrix
220    def get_q_weight_1d(self):
221        """Return 1D quadrature weights matrix of a tensor product Basis.
222
223           Returns:
224             *array: Numpy array"""
225
226        # Compute the length of the array
227        nqpts_pointer = ffi.new("CeedInt *")
228        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
229        length = nqpts_pointer[0]
230
231        # Setup the pointer's pointer
232        array_pointer = ffi.new("CeedScalar **")
233
234        # libCEED call
235        lib.CeedBasisGetQWeights(self._pointer[0], array_pointer)
236
237        # Return array created from buffer
238        # Create buffer object from returned pointer
239        buff = ffi.buffer(
240            array_pointer[0],
241            ffi.sizeof("CeedScalar") *
242            length)
243        # return read only Numpy array
244        ret = np.frombuffer(buff, dtype="float64")
245        ret.flags['WRITEABLE'] = False
246        return ret
247
248    # Get 1D quadrature points matrix
249    def get_q_ref_1d(self):
250        """Return 1D quadrature points matrix of a tensor product Basis.
251
252           Returns:
253             *array: Numpy array"""
254
255        # Compute the length of the array
256        nqpts_pointer = ffi.new("CeedInt *")
257        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
258        length = nqpts_pointer[0]
259
260        # Setup the pointer's pointer
261        array_pointer = ffi.new("CeedScalar **")
262
263        # libCEED call
264        lib.CeedBasisGetQRef(self._pointer[0], array_pointer)
265
266        # Return array created from buffer
267        # Create buffer object from returned pointer
268        buff = ffi.buffer(
269            array_pointer[0],
270            ffi.sizeof("CeedScalar") *
271            length)
272        # return read only Numpy array
273        ret = np.frombuffer(buff, dtype="float64")
274        ret.flags['WRITEABLE'] = False
275        return ret
276
277
278# ------------------------------------------------------------------------------
279
280
281class BasisTensorH1Lagrange(BasisTensorH1):
282    """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange basis
283         objects for H^1 discretizations."""
284
285    # Constructor
286    def __init__(self, ceed, dim, ncomp, P, Q, qmode):
287
288        # Setup arguments
289        self._pointer = ffi.new("CeedBasis *")
290
291        self._ceed = ceed
292
293        # libCEED call
294        err_code = lib.CeedBasisCreateTensorH1Lagrange(self._ceed._pointer[0], dim,
295                                                       ncomp, P, Q, qmode, self._pointer)
296        self._ceed._check_error(err_code)
297
298# ------------------------------------------------------------------------------
299
300
301class BasisH1(Basis):
302    """Ceed H1 Basis: finite element non tensor-product basis for H^1 discretizations."""
303
304    # Constructor
305    def __init__(self, ceed, topo, ncomp, nnodes,
306                 nqpts, interp, grad, qref, qweight):
307
308        # Setup arguments
309        self._pointer = ffi.new("CeedBasis *")
310
311        self._ceed = ceed
312
313        interp_pointer = ffi.new("CeedScalar *")
314        interp_pointer = ffi.cast(
315            "CeedScalar *",
316            interp.__array_interface__['data'][0])
317
318        grad_pointer = ffi.new("CeedScalar *")
319        grad_pointer = ffi.cast(
320            "CeedScalar *",
321            grad.__array_interface__['data'][0])
322
323        qref_pointer = ffi.new("CeedScalar *")
324        qref_pointer = ffi.cast(
325            "CeedScalar *",
326            qref.__array_interface__['data'][0])
327
328        qweight_pointer = ffi.new("CeedScalar *")
329        qweight_pointer = ffi.cast(
330            "CeedScalar *",
331            qweight.__array_interface__['data'][0])
332
333        # libCEED call
334        err_code = lib.CeedBasisCreateH1(self._ceed._pointer[0], topo, ncomp,
335                                         nnodes, nqpts, interp_pointer,
336                                         grad_pointer, qref_pointer,
337                                         qweight_pointer, self._pointer)
338
339# ------------------------------------------------------------------------------
340
341
342class BasisHdiv(Basis):
343    """Ceed Hdiv Basis: finite element non tensor-product basis for H(div) discretizations."""
344
345    # Constructor
346    def __init__(self, ceed, topo, ncomp, nnodes,
347                 nqpts, interp, grad, qref, qweight):
348
349        # Setup arguments
350        self._pointer = ffi.new("CeedBasis *")
351
352        self._ceed = ceed
353
354        interp_pointer = ffi.new("CeedScalar *")
355        interp_pointer = ffi.cast(
356            "CeedScalar *",
357            interp.__array_interface__['data'][0])
358
359        div_pointer = ffi.new("CeedScalar *")
360        div_pointer = ffi.cast(
361            "CeedScalar *",
362            grad.__array_interface__['data'][0])
363
364        qref_pointer = ffi.new("CeedScalar *")
365        qref_pointer = ffi.cast(
366            "CeedScalar *",
367            qref.__array_interface__['data'][0])
368
369        qweight_pointer = ffi.new("CeedScalar *")
370        qweight_pointer = ffi.cast(
371            "CeedScalar *",
372            qweight.__array_interface__['data'][0])
373
374        # libCEED call
375        err_code = lib.CeedBasisCreateHdiv(self._ceed._pointer[0], topo, ncomp,
376                                           nnodes, nqpts, interp_pointer,
377                                           div_pointer, qref_pointer,
378                                           qweight_pointer, self._pointer)
379
380# ------------------------------------------------------------------------------
381
382
383class BasisHcurl(Basis):
384    """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl) discretizations."""
385
386    # Constructor
387    def __init__(self, ceed, topo, ncomp, nnodes,
388                 nqpts, interp, grad, qref, qweight):
389
390        # Setup arguments
391        self._pointer = ffi.new("CeedBasis *")
392
393        self._ceed = ceed
394
395        interp_pointer = ffi.new("CeedScalar *")
396        interp_pointer = ffi.cast(
397            "CeedScalar *",
398            interp.__array_interface__['data'][0])
399
400        curl_pointer = ffi.new("CeedScalar *")
401        curl_pointer = ffi.cast(
402            "CeedScalar *",
403            grad.__array_interface__['data'][0])
404
405        qref_pointer = ffi.new("CeedScalar *")
406        qref_pointer = ffi.cast(
407            "CeedScalar *",
408            qref.__array_interface__['data'][0])
409
410        qweight_pointer = ffi.new("CeedScalar *")
411        qweight_pointer = ffi.cast(
412            "CeedScalar *",
413            qweight.__array_interface__['data'][0])
414
415        # libCEED call
416        err_code = lib.CeedBasisCreateHcurl(self._ceed._pointer[0], topo, ncomp,
417                                            nnodes, nqpts, interp_pointer,
418                                            curl_pointer, qref_pointer,
419                                            qweight_pointer, self._pointer)
420
421# ------------------------------------------------------------------------------
422
423
424class TransposeBasis():
425    """Transpose Ceed Basis: transpose of finite element tensor-product basis objects."""
426
427    # Attributes
428    _basis = None
429
430    # Constructor
431    def __init__(self, basis):
432
433        # Reference basis
434        self._basis = basis
435
436    # Representation
437    def __repr__(self):
438        return "<Transpose CeedBasis instance at " + hex(id(self)) + ">"
439
440    # Apply Transpose Basis
441    def apply(self, nelem, emode, u, v):
442        """Apply basis evaluation from quadrature points to nodes.
443
444           Args:
445             nelem: the number of elements to apply the basis evaluation to;
446                      the backend will specify the ordering in a
447                      Blocked ElemRestriction
448             **emode: basis evaluation mode
449             u: input vector
450             v: output vector"""
451
452        # libCEED call
453        self._basis.apply(nelem, emode, u, v, tmode=TRANSPOSE)
454
455# ------------------------------------------------------------------------------
456