xref: /libCEED/python/ceed_elemrestriction.py (revision 187168c7adf497522cc3b03930f59faf3536934a)
1c8e769f6Sjeremylt# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2c8e769f6Sjeremylt# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3c8e769f6Sjeremylt# reserved. See files LICENSE and NOTICE for details.
4c8e769f6Sjeremylt#
5c8e769f6Sjeremylt# This file is part of CEED, a collection of benchmarks, miniapps, software
6c8e769f6Sjeremylt# libraries and APIs for efficient high-order finite element and spectral
7c8e769f6Sjeremylt# element discretizations for exascale applications. For more information and
8c8e769f6Sjeremylt# source code availability see http://github.com/ceed.
9c8e769f6Sjeremylt#
10c8e769f6Sjeremylt# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11c8e769f6Sjeremylt# a collaborative effort of two U.S. Department of Energy organizations (Office
12c8e769f6Sjeremylt# of Science and the National Nuclear Security Administration) responsible for
13c8e769f6Sjeremylt# the planning and preparation of a capable exascale ecosystem, including
14c8e769f6Sjeremylt# software, applications, hardware, advanced system engineering and early
15c8e769f6Sjeremylt# testbed platforms, in support of the nation's exascale computing imperative.
16c8e769f6Sjeremylt
17c8e769f6Sjeremyltfrom _ceed_cffi import ffi, lib
18c8e769f6Sjeremyltimport tempfile
19c8e769f6Sjeremyltimport numpy as np
20c8e769f6Sjeremyltfrom abc import ABC
21*187168c7SJeremy L Thompsonfrom .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, MEM_HOST, USE_POINTER, COPY_VALUES, TRANSPOSE, NOTRANSPOSE
22c8e769f6Sjeremyltfrom .ceed_vector import _VectorWrap
23c8e769f6Sjeremylt
24c8e769f6Sjeremylt# ------------------------------------------------------------------------------
257a7b0fa3SJed Brown
267a7b0fa3SJed Brown
27c8e769f6Sjeremyltclass _ElemRestrictionBase(ABC):
28c8e769f6Sjeremylt    """Ceed ElemRestriction: restriction from local vectors to elements."""
29c8e769f6Sjeremylt
30c8e769f6Sjeremylt    # Attributes
31*187168c7SJeremy L Thompson    _ceed = None
32c8e769f6Sjeremylt    _pointer = ffi.NULL
33*187168c7SJeremy L Thompson    _array_reference = None
34c8e769f6Sjeremylt
35c8e769f6Sjeremylt    # Destructor
36c8e769f6Sjeremylt    def __del__(self):
37c8e769f6Sjeremylt        # libCEED call
38477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionDestroy(self._pointer)
39477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
40c8e769f6Sjeremylt
41c8e769f6Sjeremylt    # Representation
42c8e769f6Sjeremylt    def __repr__(self):
43c8e769f6Sjeremylt        return "<CeedElemRestriction instance at " + hex(id(self)) + ">"
44c8e769f6Sjeremylt
45c8e769f6Sjeremylt    # String conversion for print() to stdout
46c8e769f6Sjeremylt    def __str__(self):
47c8e769f6Sjeremylt        """View an ElemRestriction via print()."""
48c8e769f6Sjeremylt
49c8e769f6Sjeremylt        # libCEED call
50c8e769f6Sjeremylt        with tempfile.NamedTemporaryFile() as key_file:
51c8e769f6Sjeremylt            with open(key_file.name, 'r+') as stream_file:
52c8e769f6Sjeremylt                stream = ffi.cast("FILE *", stream_file)
53c8e769f6Sjeremylt
54477729cfSJeremy L Thompson                err_code = lib.CeedElemRestrictionView(self._pointer[0], stream)
55477729cfSJeremy L Thompson                self._ceed._check_error(err_code)
56c8e769f6Sjeremylt
57c8e769f6Sjeremylt                stream_file.seek(0)
58c8e769f6Sjeremylt                out_string = stream_file.read()
59c8e769f6Sjeremylt
60c8e769f6Sjeremylt        return out_string
61c8e769f6Sjeremylt
62c8e769f6Sjeremylt    # Apply CeedElemRestriction
63a8d32208Sjeremylt    def apply(self, u, v, tmode=NOTRANSPOSE, request=REQUEST_IMMEDIATE):
64c8e769f6Sjeremylt        """Restrict a local vector to an element vector or apply its transpose.
65c8e769f6Sjeremylt
66c8e769f6Sjeremylt           Args:
67c8e769f6Sjeremylt             u: input vector
68c8e769f6Sjeremylt             v: output vector
69c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
70c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
71c8e769f6Sjeremylt
72c8e769f6Sjeremylt        # libCEED call
73477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionApply(self._pointer[0], tmode, u._pointer[0],
74a8d32208Sjeremylt                                                v._pointer[0], request)
75477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
76c8e769f6Sjeremylt
77c8e769f6Sjeremylt    # Transpose an ElemRestriction
78c8e769f6Sjeremylt    @property
79c8e769f6Sjeremylt    def T(self):
80c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
81c8e769f6Sjeremylt
82c8e769f6Sjeremylt        return TransposeElemRestriction(self)
83c8e769f6Sjeremylt
84c8e769f6Sjeremylt    # Transpose an ElemRestriction
85c8e769f6Sjeremylt    @property
86c8e769f6Sjeremylt    def transpose(self):
87c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
88c8e769f6Sjeremylt
89c8e769f6Sjeremylt        return TransposeElemRestriction(self)
90c8e769f6Sjeremylt
91c8e769f6Sjeremylt    # Create restriction vectors
92c8e769f6Sjeremylt    def create_vector(self, createLvec=True, createEvec=True):
93c8e769f6Sjeremylt        """Create Vectors associated with an ElemRestriction.
94c8e769f6Sjeremylt
95c8e769f6Sjeremylt           Args:
96c8e769f6Sjeremylt             **createLvec: flag to create local vector, default True
97c8e769f6Sjeremylt             **createEvec: flag to create element vector, default True
98c8e769f6Sjeremylt
99c8e769f6Sjeremylt           Returns:
100c8e769f6Sjeremylt             [lvec, evec]: local vector and element vector, or None if flag set to false"""
101c8e769f6Sjeremylt
102c8e769f6Sjeremylt        # Vector pointers
103c8e769f6Sjeremylt        lvecPointer = ffi.new("CeedVector *") if createLvec else ffi.NULL
104c8e769f6Sjeremylt        evecPointer = ffi.new("CeedVector *") if createEvec else ffi.NULL
105c8e769f6Sjeremylt
106c8e769f6Sjeremylt        # libCEED call
107477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateVector(self._pointer[0], lvecPointer,
108c8e769f6Sjeremylt                                                       evecPointer)
109477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
110c8e769f6Sjeremylt
111c8e769f6Sjeremylt        # Return vectors
1127a7b0fa3SJed Brown        lvec = _VectorWrap(
113477729cfSJeremy L Thompson            self._ceed, lvecPointer) if createLvec else None
1147a7b0fa3SJed Brown        evec = _VectorWrap(
115477729cfSJeremy L Thompson            self._ceed, evecPointer) if createEvec else None
116c8e769f6Sjeremylt
117c8e769f6Sjeremylt        # Return
118c8e769f6Sjeremylt        return [lvec, evec]
119c8e769f6Sjeremylt
120c8e769f6Sjeremylt    # Get ElemRestriction multiplicity
121a8d32208Sjeremylt    def get_multiplicity(self):
122c8e769f6Sjeremylt        """Get the multiplicity of nodes in an ElemRestriction.
123c8e769f6Sjeremylt
124c8e769f6Sjeremylt           Returns:
125c8e769f6Sjeremylt             mult: local vector containing multiplicity of nodes in ElemRestriction"""
126c8e769f6Sjeremylt
127c8e769f6Sjeremylt        # Create mult vector
128c8e769f6Sjeremylt        [mult, evec] = self.create_vector(createEvec=False)
129c8e769f6Sjeremylt        mult.set_value(0)
130c8e769f6Sjeremylt
131c8e769f6Sjeremylt        # libCEED call
132477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionGetMultiplicity(
1337a7b0fa3SJed Brown            self._pointer[0], mult._pointer[0])
134477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
135c8e769f6Sjeremylt
136c8e769f6Sjeremylt        # Return
137c8e769f6Sjeremylt        return mult
138c8e769f6Sjeremylt
139c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1407a7b0fa3SJed Brown
1417a7b0fa3SJed Brown
142c8e769f6Sjeremyltclass ElemRestriction(_ElemRestrictionBase):
143c8e769f6Sjeremylt    """Ceed ElemRestriction: restriction from local vectors to elements."""
144c8e769f6Sjeremylt
145c8e769f6Sjeremylt    # Constructor
146d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets,
147d979a051Sjeremylt                 memtype=MEM_HOST, cmode=COPY_VALUES):
148c8e769f6Sjeremylt        # CeedVector object
149c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
150c8e769f6Sjeremylt
151c8e769f6Sjeremylt        # Reference to Ceed
152c8e769f6Sjeremylt        self._ceed = ceed
153c8e769f6Sjeremylt
154*187168c7SJeremy L Thompson        # Store array reference if needed
155*187168c7SJeremy L Thompson        if cmode == USE_POINTER:
156*187168c7SJeremy L Thompson            self._array_reference = offsets
157*187168c7SJeremy L Thompson        else:
158*187168c7SJeremy L Thompson            self._array_reference = None
159*187168c7SJeremy L Thompson
160c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
161d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
162d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
163d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
164c8e769f6Sjeremylt
165c8e769f6Sjeremylt        # libCEED call
166477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreate(self._ceed._pointer[0], nelem,
167477729cfSJeremy L Thompson                                                 elemsize, ncomp, compstride,
168477729cfSJeremy L Thompson                                                 lsize, memtype, cmode,
169d979a051Sjeremylt                                                 offsets_pointer, self._pointer)
170477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
171c8e769f6Sjeremylt
172c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1737a7b0fa3SJed Brown
1747a7b0fa3SJed Brown
17569a53589Sjeremyltclass StridedElemRestriction(_ElemRestrictionBase):
17669a53589Sjeremylt    """Ceed Strided ElemRestriction: strided restriction from local vectors to elements."""
177c8e769f6Sjeremylt
178c8e769f6Sjeremylt    # Constructor
179d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, lsize, strides):
180c8e769f6Sjeremylt        # CeedVector object
181c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
182c8e769f6Sjeremylt
183c8e769f6Sjeremylt        # Reference to Ceed
184c8e769f6Sjeremylt        self._ceed = ceed
185c8e769f6Sjeremylt
18669a53589Sjeremylt        # Setup the numpy array for the libCEED call
18769a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
18869a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
18969a53589Sjeremylt                                   strides.__array_interface__['data'][0])
19069a53589Sjeremylt
191c8e769f6Sjeremylt        # libCEED call
192477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0],
193477729cfSJeremy L Thompson                                                        nelem, elemsize, ncomp,
194477729cfSJeremy L Thompson                                                        lsize, strides_pointer,
195477729cfSJeremy L Thompson                                                        self._pointer)
196477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
197c8e769f6Sjeremylt
198c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1997a7b0fa3SJed Brown
2007a7b0fa3SJed Brown
201c8e769f6Sjeremyltclass BlockedElemRestriction(_ElemRestrictionBase):
202c8e769f6Sjeremylt    """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements."""
203c8e769f6Sjeremylt
204c8e769f6Sjeremylt    # Constructor
205d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize,
206d979a051Sjeremylt                 offsets, memtype=MEM_HOST, cmode=COPY_VALUES):
207c8e769f6Sjeremylt        # CeedVector object
208c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
209c8e769f6Sjeremylt
210c8e769f6Sjeremylt        # Reference to Ceed
211c8e769f6Sjeremylt        self._ceed = ceed
212c8e769f6Sjeremylt
213c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
214d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
215d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
216d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
217c8e769f6Sjeremylt
218c8e769f6Sjeremylt        # libCEED call
219477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], nelem,
220d979a051Sjeremylt                                                        elemsize, blksize, ncomp,
221d979a051Sjeremylt                                                        compstride, lsize, memtype, cmode,
222d979a051Sjeremylt                                                        offsets_pointer, self._pointer)
223477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
224c8e769f6Sjeremylt
225c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
226c8e769f6Sjeremylt    @property
227c8e769f6Sjeremylt    def T(self):
228c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
229c8e769f6Sjeremylt
230c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
231c8e769f6Sjeremylt
232c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
233c8e769f6Sjeremylt    @property
234c8e769f6Sjeremylt    def transpose(self):
235c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
236c8e769f6Sjeremylt
237c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
238c8e769f6Sjeremylt
239c8e769f6Sjeremylt    # Apply CeedElemRestriction to single block
240a8d32208Sjeremylt    def apply_block(self, block, u, v, tmode=NOTRANSPOSE,
241c8e769f6Sjeremylt                    request=REQUEST_IMMEDIATE):
242c8e769f6Sjeremylt        """Restrict a local vector to a block of an element vector or apply its transpose.
243c8e769f6Sjeremylt
244c8e769f6Sjeremylt           Args:
245c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
246c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
247c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
248c8e769f6Sjeremylt             u: input vector
249c8e769f6Sjeremylt             v: output vector
250c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
251c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
252c8e769f6Sjeremylt
253c8e769f6Sjeremylt        # libCEED call
254477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode,
255477729cfSJeremy L Thompson                                                     u._pointer[0], v._pointer[0],
256477729cfSJeremy L Thompson                                                     request)
257477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
258c8e769f6Sjeremylt
259c8e769f6Sjeremylt# ------------------------------------------------------------------------------
2607a7b0fa3SJed Brown
2617a7b0fa3SJed Brown
26269a53589Sjeremyltclass BlockedStridedElemRestriction(BlockedElemRestriction):
26369a53589Sjeremylt    """Ceed Blocked Strided ElemRestriction: strided restriction from local vectors to elements."""
26469a53589Sjeremylt
26569a53589Sjeremylt    # Constructor
266d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, lsize, strides):
26769a53589Sjeremylt        # CeedVector object
26869a53589Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
26969a53589Sjeremylt
27069a53589Sjeremylt        # Reference to Ceed
27169a53589Sjeremylt        self._ceed = ceed
27269a53589Sjeremylt
27369a53589Sjeremylt        # Setup the numpy array for the libCEED call
27469a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
27569a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
27669a53589Sjeremylt                                   strides.__array_interface__['data'][0])
27769a53589Sjeremylt
27869a53589Sjeremylt        # libCEED call
279477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem,
280d979a051Sjeremylt                                                               elemsize, blksize, ncomp,
281d979a051Sjeremylt                                                               lsize, strides_pointer,
28269a53589Sjeremylt                                                               self._pointer)
283477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
28469a53589Sjeremylt
28569a53589Sjeremylt# ------------------------------------------------------------------------------
2867a7b0fa3SJed Brown
2877a7b0fa3SJed Brown
288c8e769f6Sjeremyltclass TransposeElemRestriction():
289c8e769f6Sjeremylt    """Ceed ElemRestriction: transpose restriction from elements to local vectors."""
290c8e769f6Sjeremylt
291c8e769f6Sjeremylt    # Attributes
292c8e769f6Sjeremylt    _elemrestriction = None
293c8e769f6Sjeremylt
294c8e769f6Sjeremylt    # Constructor
295c8e769f6Sjeremylt    def __init__(self, elemrestriction):
296c8e769f6Sjeremylt
297c8e769f6Sjeremylt        # Reference elemrestriction
298c8e769f6Sjeremylt        self._elemrestriction = elemrestriction
299c8e769f6Sjeremylt
300c8e769f6Sjeremylt    # Representation
301c8e769f6Sjeremylt    def __repr__(self):
3027a7b0fa3SJed Brown        return "<Transpose CeedElemRestriction instance at " + \
3037a7b0fa3SJed Brown            hex(id(self)) + ">"
304c8e769f6Sjeremylt
305c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
3067a7b0fa3SJed Brown
307a8d32208Sjeremylt    def apply(self, u, v, request=REQUEST_IMMEDIATE):
308c8e769f6Sjeremylt        """Restrict an element vector to a local vector.
309c8e769f6Sjeremylt
310c8e769f6Sjeremylt           Args:
311c8e769f6Sjeremylt             u: input vector
312c8e769f6Sjeremylt             v: output vector
313c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
314c8e769f6Sjeremylt
315c8e769f6Sjeremylt        # libCEED call
316a8d32208Sjeremylt        self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE)
317c8e769f6Sjeremylt
318c8e769f6Sjeremylt# ------------------------------------------------------------------------------
3197a7b0fa3SJed Brown
3207a7b0fa3SJed Brown
321c8e769f6Sjeremyltclass TransposeBlockedElemRestriction(TransposeElemRestriction):
322c8e769f6Sjeremylt    """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements
323c8e769f6Sjeremylt         to local vectors."""
324c8e769f6Sjeremylt
325c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
326a8d32208Sjeremylt    def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE):
327c8e769f6Sjeremylt        """Restrict a block of an element vector to a local vector.
328c8e769f6Sjeremylt
329c8e769f6Sjeremylt           Args:
330c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
331c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
332c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
333c8e769f6Sjeremylt             u: input vector
334c8e769f6Sjeremylt             v: output vector
335c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
336c8e769f6Sjeremylt
337c8e769f6Sjeremylt        # libCEED call
338a8d32208Sjeremylt        self._elemrestriction.apply_block(block, u, v, request=request,
339c8e769f6Sjeremylt                                          tmode=TRANSPOSE)
340c8e769f6Sjeremylt
341c8e769f6Sjeremylt# ------------------------------------------------------------------------------
342