xref: /libCEED/python/ceed_basis.py (revision 20e46440c131be8f30f4ad01c95e1bd89e184efa) 
1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3# reserved. See files LICENSE and NOTICE for details.
4#
5# This file is part of CEED, a collection of benchmarks, miniapps, software
6# libraries and APIs for efficient high-order finite element and spectral
7# element discretizations for exascale applications. For more information and
8# source code availability see http://github.com/ceed.
9#
10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11# a collaborative effort of two U.S. Department of Energy organizations (Office
12# of Science and the National Nuclear Security Administration) responsible for
13# the planning and preparation of a capable exascale ecosystem, including
14# software, applications, hardware, advanced system engineering and early
15# testbed platforms, in support of the nation's exascale computing imperative.
16
17from _ceed_cffi import ffi, lib
18import tempfile
19import numpy as np
20from abc import ABC
21from .ceed_constants import TRANSPOSE, NOTRANSPOSE, scalar_types
22
23# ------------------------------------------------------------------------------
24
25
26class Basis(ABC):
27    """Ceed Basis: finite element basis objects."""
28
29    # Representation
30    def __repr__(self):
31        return "<CeedBasis instance at " + hex(id(self)) + ">"
32
33    # String conversion for print() to stdout
34    def __str__(self):
35        """View a Basis via print()."""
36
37        # libCEED call
38        with tempfile.NamedTemporaryFile() as key_file:
39            with open(key_file.name, 'r+') as stream_file:
40                stream = ffi.cast("FILE *", stream_file)
41
42                err_code = lib.CeedBasisView(self._pointer[0], stream)
43                self._ceed._check_error(err_code)
44
45                stream_file.seek(0)
46                out_string = stream_file.read()
47
48        return out_string
49
50    # Apply Basis
51    def apply(self, nelem, emode, u, v, tmode=NOTRANSPOSE):
52        """Apply basis evaluation from nodes to quadrature points or vice versa.
53
54           Args:
55             nelem: the number of elements to apply the basis evaluation to;
56                      the backend will specify the ordering in a
57                      BlockedElemRestriction
58             emode: basis evaluation mode
59             u: input vector
60             v: output vector
61             **tmode: CEED_NOTRANSPOSE to evaluate from nodes to quadrature
62                        points, CEED_TRANSPOSE to apply the transpose, mapping
63                        from quadrature points to nodes; default CEED_NOTRANSPOSE"""
64
65        # libCEED call
66        err_code = lib.CeedBasisApply(self._pointer[0], nelem, tmode, emode,
67                                      u._pointer[0], v._pointer[0])
68        self._ceed._check_error(err_code)
69
70    # Transpose a Basis
71    @property
72    def T(self):
73        """Transpose a Basis."""
74
75        return TransposeBasis(self)
76
77    # Transpose a Basis
78    @property
79    def transpose(self):
80        """Transpose a Basis."""
81
82        return TransposeBasis(self)
83
84    # Get number of nodes
85    def get_num_nodes(self):
86        """Get total number of nodes (in dim dimensions) of a Basis.
87
88           Returns:
89             num_nodes: total number of nodes"""
90
91        # Setup argument
92        p_pointer = ffi.new("CeedInt *")
93
94        # libCEED call
95        err_code = lib.CeedBasisGetNumNodes(self._pointer[0], p_pointer)
96        self._ceed._check_error(err_code)
97
98        return p_pointer[0]
99
100    # Get number of quadrature points
101    def get_num_quadrature_points(self):
102        """Get total number of quadrature points (in dim dimensions) of a Basis.
103
104           Returns:
105             num_qpts: total number of quadrature points"""
106
107        # Setup argument
108        q_pointer = ffi.new("CeedInt *")
109
110        # libCEED call
111        err_code = lib.CeedBasisGetNumQuadraturePoints(
112            self._pointer[0], q_pointer)
113        self._ceed._check_error(err_code)
114
115        return q_pointer[0]
116
117    # Destructor
118    def __del__(self):
119        # libCEED call
120        err_code = lib.CeedBasisDestroy(self._pointer)
121        self._ceed._check_error(err_code)
122
123# ------------------------------------------------------------------------------
124
125
126class BasisTensorH1(Basis):
127    """Ceed Tensor H1 Basis: finite element tensor-product basis objects for
128         H^1 discretizations."""
129
130    # Constructor
131    def __init__(self, ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d,
132                 qref1d, qweight1d):
133
134        # Setup arguments
135        self._pointer = ffi.new("CeedBasis *")
136
137        self._ceed = ceed
138
139        interp1d_pointer = ffi.new("CeedScalar *")
140        interp1d_pointer = ffi.cast(
141            "CeedScalar *",
142            interp1d.__array_interface__['data'][0])
143
144        grad1d_pointer = ffi.new("CeedScalar *")
145        grad1d_pointer = ffi.cast(
146            "CeedScalar *",
147            grad1d.__array_interface__['data'][0])
148
149        qref1d_pointer = ffi.new("CeedScalar *")
150        qref1d_pointer = ffi.cast(
151            "CeedScalar *",
152            qref1d.__array_interface__['data'][0])
153
154        qweight1d_pointer = ffi.new("CeedScalar *")
155        qweight1d_pointer = ffi.cast(
156            "CeedScalar *",
157            qweight1d.__array_interface__['data'][0])
158
159        # libCEED call
160        err_code = lib.CeedBasisCreateTensorH1(self._ceed._pointer[0], dim, ncomp,
161                                               P1d, Q1d, interp1d_pointer,
162                                               grad1d_pointer, qref1d_pointer,
163                                               qweight1d_pointer, self._pointer)
164        self._ceed._check_error(err_code)
165
166    # Get 1D interpolation matrix
167    def get_interp_1d(self):
168        """Return 1D interpolation matrix of a tensor product Basis.
169
170           Returns:
171             *array: Numpy array"""
172
173        # Compute the length of the array
174        nnodes_pointer = ffi.new("CeedInt *")
175        lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer)
176        nqpts_pointer = ffi.new("CeedInt *")
177        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
178        length = nnodes_pointer[0] * nqpts_pointer[0]
179
180        # Setup the pointer's pointer
181        array_pointer = ffi.new("CeedScalar **")
182
183        # libCEED call
184        lib.CeedBasisGetInterp1D(self._pointer[0], array_pointer)
185
186        # Return array created from buffer
187        # Create buffer object from returned pointer
188        buff = ffi.buffer(
189            array_pointer[0],
190            ffi.sizeof("CeedScalar") *
191            length)
192        # return read only Numpy array
193        ret = np.frombuffer(buff, dtype=scalar_types[lib.CEED_SCALAR_TYPE])
194        ret.flags['WRITEABLE'] = False
195        return ret
196
197    # Get 1D gradient matrix
198    def get_grad_1d(self):
199        """Return 1D gradient matrix of a tensor product Basis.
200
201           Returns:
202             *array: Numpy array"""
203
204        # Compute the length of the array
205        nnodes_pointer = ffi.new("CeedInt *")
206        lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer)
207        nqpts_pointer = ffi.new("CeedInt *")
208        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
209        length = nnodes_pointer[0] * nqpts_pointer[0]
210
211        # Setup the pointer's pointer
212        array_pointer = ffi.new("CeedScalar **")
213
214        # libCEED call
215        lib.CeedBasisGetGrad1D(self._pointer[0], array_pointer)
216
217        # Return array created from buffer
218        # Create buffer object from returned pointer
219        buff = ffi.buffer(
220            array_pointer[0],
221            ffi.sizeof("CeedScalar") *
222            length)
223        # return read only Numpy array
224        ret = np.frombuffer(buff, dtype="float64")
225        ret.flags['WRITEABLE'] = False
226        return ret
227
228    # Get 1D quadrature weights matrix
229    def get_q_weight_1d(self):
230        """Return 1D quadrature weights matrix of a tensor product Basis.
231
232           Returns:
233             *array: Numpy array"""
234
235        # Compute the length of the array
236        nqpts_pointer = ffi.new("CeedInt *")
237        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
238        length = nqpts_pointer[0]
239
240        # Setup the pointer's pointer
241        array_pointer = ffi.new("CeedScalar **")
242
243        # libCEED call
244        lib.CeedBasisGetQWeights(self._pointer[0], array_pointer)
245
246        # Return array created from buffer
247        # Create buffer object from returned pointer
248        buff = ffi.buffer(
249            array_pointer[0],
250            ffi.sizeof("CeedScalar") *
251            length)
252        # return read only Numpy array
253        ret = np.frombuffer(buff, dtype="float64")
254        ret.flags['WRITEABLE'] = False
255        return ret
256
257    # Get 1D quadrature points matrix
258    def get_q_ref_1d(self):
259        """Return 1D quadrature points matrix of a tensor product Basis.
260
261           Returns:
262             *array: Numpy array"""
263
264        # Compute the length of the array
265        nqpts_pointer = ffi.new("CeedInt *")
266        lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer)
267        length = nqpts_pointer[0]
268
269        # Setup the pointer's pointer
270        array_pointer = ffi.new("CeedScalar **")
271
272        # libCEED call
273        lib.CeedBasisGetQRef(self._pointer[0], array_pointer)
274
275        # Return array created from buffer
276        # Create buffer object from returned pointer
277        buff = ffi.buffer(
278            array_pointer[0],
279            ffi.sizeof("CeedScalar") *
280            length)
281        # return read only Numpy array
282        ret = np.frombuffer(buff, dtype="float64")
283        ret.flags['WRITEABLE'] = False
284        return ret
285
286
287# ------------------------------------------------------------------------------
288
289
290class BasisTensorH1Lagrange(BasisTensorH1):
291    """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange basis
292         objects for H^1 discretizations."""
293
294    # Constructor
295    def __init__(self, ceed, dim, ncomp, P, Q, qmode):
296
297        # Setup arguments
298        self._pointer = ffi.new("CeedBasis *")
299
300        self._ceed = ceed
301
302        # libCEED call
303        err_code = lib.CeedBasisCreateTensorH1Lagrange(self._ceed._pointer[0], dim,
304                                                       ncomp, P, Q, qmode, self._pointer)
305        self._ceed._check_error(err_code)
306
307# ------------------------------------------------------------------------------
308
309
310class BasisH1(Basis):
311    """Ceed H1 Basis: finite element non tensor-product basis for H^1 discretizations."""
312
313    # Constructor
314    def __init__(self, ceed, topo, ncomp, nnodes,
315                 nqpts, interp, grad, qref, qweight):
316
317        # Setup arguments
318        self._pointer = ffi.new("CeedBasis *")
319
320        self._ceed = ceed
321
322        interp_pointer = ffi.new("CeedScalar *")
323        interp_pointer = ffi.cast(
324            "CeedScalar *",
325            interp.__array_interface__['data'][0])
326
327        grad_pointer = ffi.new("CeedScalar *")
328        grad_pointer = ffi.cast(
329            "CeedScalar *",
330            grad.__array_interface__['data'][0])
331
332        qref_pointer = ffi.new("CeedScalar *")
333        qref_pointer = ffi.cast(
334            "CeedScalar *",
335            qref.__array_interface__['data'][0])
336
337        qweight_pointer = ffi.new("CeedScalar *")
338        qweight_pointer = ffi.cast(
339            "CeedScalar *",
340            qweight.__array_interface__['data'][0])
341
342        # libCEED call
343        err_code = lib.CeedBasisCreateH1(self._ceed._pointer[0], topo, ncomp,
344                                         nnodes, nqpts, interp_pointer,
345                                         grad_pointer, qref_pointer,
346                                         qweight_pointer, self._pointer)
347
348# ------------------------------------------------------------------------------
349
350
351class TransposeBasis():
352    """Transpose Ceed Basis: transpose of finite element tensor-product basis objects."""
353
354    # Attributes
355    _basis = None
356
357    # Constructor
358    def __init__(self, basis):
359
360        # Reference basis
361        self._basis = basis
362
363    # Representation
364    def __repr__(self):
365        return "<Transpose CeedBasis instance at " + hex(id(self)) + ">"
366
367    # Apply Transpose Basis
368    def apply(self, nelem, emode, u, v):
369        """Apply basis evaluation from quadrature points to nodes.
370
371           Args:
372             nelem: the number of elements to apply the basis evaluation to;
373                      the backend will specify the ordering in a
374                      Blocked ElemRestriction
375             **emode: basis evaluation mode
376             u: input vector
377             v: output vector"""
378
379        # libCEED call
380        self._basis.apply(nelem, emode, u, v, tmode=TRANSPOSE)
381
382# ------------------------------------------------------------------------------
383