xref: /libCEED/python/ceed_basis.py (revision ce18bed930e8f3bfebcf709a18844aba97fe4630) 
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 TransposeBasis():
343    """Transpose Ceed Basis: transpose of finite element tensor-product basis objects."""
344
345    # Attributes
346    _basis = None
347
348    # Constructor
349    def __init__(self, basis):
350
351        # Reference basis
352        self._basis = basis
353
354    # Representation
355    def __repr__(self):
356        return "<Transpose CeedBasis instance at " + hex(id(self)) + ">"
357
358    # Apply Transpose Basis
359    def apply(self, nelem, emode, u, v):
360        """Apply basis evaluation from quadrature points to nodes.
361
362           Args:
363             nelem: the number of elements to apply the basis evaluation to;
364                      the backend will specify the ordering in a
365                      Blocked ElemRestriction
366             **emode: basis evaluation mode
367             u: input vector
368             v: output vector"""
369
370        # libCEED call
371        self._basis.apply(nelem, emode, u, v, tmode=TRANSPOSE)
372
373# ------------------------------------------------------------------------------
374