xref: /libCEED/python/ceed_elemrestriction.py (revision 477729cfdbe3383ce841609affbc919aa37fc70d)
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
21d979a051Sjeremyltfrom .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, MEM_HOST, 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
31c8e769f6Sjeremylt    _ceed = ffi.NULL
32c8e769f6Sjeremylt    _pointer = ffi.NULL
33c8e769f6Sjeremylt
34c8e769f6Sjeremylt    # Destructor
35c8e769f6Sjeremylt    def __del__(self):
36c8e769f6Sjeremylt        # libCEED call
37*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionDestroy(self._pointer)
38*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
39c8e769f6Sjeremylt
40c8e769f6Sjeremylt    # Representation
41c8e769f6Sjeremylt    def __repr__(self):
42c8e769f6Sjeremylt        return "<CeedElemRestriction instance at " + hex(id(self)) + ">"
43c8e769f6Sjeremylt
44c8e769f6Sjeremylt    # String conversion for print() to stdout
45c8e769f6Sjeremylt    def __str__(self):
46c8e769f6Sjeremylt        """View an ElemRestriction via print()."""
47c8e769f6Sjeremylt
48c8e769f6Sjeremylt        # libCEED call
49c8e769f6Sjeremylt        with tempfile.NamedTemporaryFile() as key_file:
50c8e769f6Sjeremylt            with open(key_file.name, 'r+') as stream_file:
51c8e769f6Sjeremylt                stream = ffi.cast("FILE *", stream_file)
52c8e769f6Sjeremylt
53*477729cfSJeremy L Thompson                err_code = lib.CeedElemRestrictionView(self._pointer[0], stream)
54*477729cfSJeremy L Thompson                self._ceed._check_error(err_code)
55c8e769f6Sjeremylt
56c8e769f6Sjeremylt                stream_file.seek(0)
57c8e769f6Sjeremylt                out_string = stream_file.read()
58c8e769f6Sjeremylt
59c8e769f6Sjeremylt        return out_string
60c8e769f6Sjeremylt
61c8e769f6Sjeremylt    # Apply CeedElemRestriction
62a8d32208Sjeremylt    def apply(self, u, v, tmode=NOTRANSPOSE, request=REQUEST_IMMEDIATE):
63c8e769f6Sjeremylt        """Restrict a local vector to an element vector or apply its transpose.
64c8e769f6Sjeremylt
65c8e769f6Sjeremylt           Args:
66c8e769f6Sjeremylt             u: input vector
67c8e769f6Sjeremylt             v: output vector
68c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
69c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
70c8e769f6Sjeremylt
71c8e769f6Sjeremylt        # libCEED call
72*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionApply(self._pointer[0], tmode, u._pointer[0],
73a8d32208Sjeremylt                                                v._pointer[0], request)
74*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
75c8e769f6Sjeremylt
76c8e769f6Sjeremylt    # Transpose an ElemRestriction
77c8e769f6Sjeremylt    @property
78c8e769f6Sjeremylt    def T(self):
79c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
80c8e769f6Sjeremylt
81c8e769f6Sjeremylt        return TransposeElemRestriction(self)
82c8e769f6Sjeremylt
83c8e769f6Sjeremylt    # Transpose an ElemRestriction
84c8e769f6Sjeremylt    @property
85c8e769f6Sjeremylt    def transpose(self):
86c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
87c8e769f6Sjeremylt
88c8e769f6Sjeremylt        return TransposeElemRestriction(self)
89c8e769f6Sjeremylt
90c8e769f6Sjeremylt    # Create restriction vectors
91c8e769f6Sjeremylt    def create_vector(self, createLvec=True, createEvec=True):
92c8e769f6Sjeremylt        """Create Vectors associated with an ElemRestriction.
93c8e769f6Sjeremylt
94c8e769f6Sjeremylt           Args:
95c8e769f6Sjeremylt             **createLvec: flag to create local vector, default True
96c8e769f6Sjeremylt             **createEvec: flag to create element vector, default True
97c8e769f6Sjeremylt
98c8e769f6Sjeremylt           Returns:
99c8e769f6Sjeremylt             [lvec, evec]: local vector and element vector, or None if flag set to false"""
100c8e769f6Sjeremylt
101c8e769f6Sjeremylt        # Vector pointers
102c8e769f6Sjeremylt        lvecPointer = ffi.new("CeedVector *") if createLvec else ffi.NULL
103c8e769f6Sjeremylt        evecPointer = ffi.new("CeedVector *") if createEvec else ffi.NULL
104c8e769f6Sjeremylt
105c8e769f6Sjeremylt        # libCEED call
106*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateVector(self._pointer[0], lvecPointer,
107c8e769f6Sjeremylt                                                       evecPointer)
108*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
109c8e769f6Sjeremylt
110c8e769f6Sjeremylt        # Return vectors
1117a7b0fa3SJed Brown        lvec = _VectorWrap(
112*477729cfSJeremy L Thompson            self._ceed, lvecPointer) if createLvec else None
1137a7b0fa3SJed Brown        evec = _VectorWrap(
114*477729cfSJeremy L Thompson            self._ceed, evecPointer) if createEvec else None
115c8e769f6Sjeremylt
116c8e769f6Sjeremylt        # Return
117c8e769f6Sjeremylt        return [lvec, evec]
118c8e769f6Sjeremylt
119c8e769f6Sjeremylt    # Get ElemRestriction multiplicity
120a8d32208Sjeremylt    def get_multiplicity(self):
121c8e769f6Sjeremylt        """Get the multiplicity of nodes in an ElemRestriction.
122c8e769f6Sjeremylt
123c8e769f6Sjeremylt           Returns:
124c8e769f6Sjeremylt             mult: local vector containing multiplicity of nodes in ElemRestriction"""
125c8e769f6Sjeremylt
126c8e769f6Sjeremylt        # Create mult vector
127c8e769f6Sjeremylt        [mult, evec] = self.create_vector(createEvec=False)
128c8e769f6Sjeremylt        mult.set_value(0)
129c8e769f6Sjeremylt
130c8e769f6Sjeremylt        # libCEED call
131*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionGetMultiplicity(
1327a7b0fa3SJed Brown            self._pointer[0], mult._pointer[0])
133*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
134c8e769f6Sjeremylt
135c8e769f6Sjeremylt        # Return
136c8e769f6Sjeremylt        return mult
137c8e769f6Sjeremylt
138c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1397a7b0fa3SJed Brown
1407a7b0fa3SJed Brown
141c8e769f6Sjeremyltclass ElemRestriction(_ElemRestrictionBase):
142c8e769f6Sjeremylt    """Ceed ElemRestriction: restriction from local vectors to elements."""
143c8e769f6Sjeremylt
144c8e769f6Sjeremylt    # Constructor
145d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets,
146d979a051Sjeremylt                 memtype=MEM_HOST, cmode=COPY_VALUES):
147c8e769f6Sjeremylt        # CeedVector object
148c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
149c8e769f6Sjeremylt
150c8e769f6Sjeremylt        # Reference to Ceed
151c8e769f6Sjeremylt        self._ceed = ceed
152c8e769f6Sjeremylt
153c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
154d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
155d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
156d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
157c8e769f6Sjeremylt
158c8e769f6Sjeremylt        # libCEED call
159*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreate(self._ceed._pointer[0], nelem,
160*477729cfSJeremy L Thompson                                                 elemsize, ncomp, compstride,
161*477729cfSJeremy L Thompson                                                 lsize, memtype, cmode,
162d979a051Sjeremylt                                                 offsets_pointer, self._pointer)
163*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
164c8e769f6Sjeremylt
165c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1667a7b0fa3SJed Brown
1677a7b0fa3SJed Brown
16869a53589Sjeremyltclass StridedElemRestriction(_ElemRestrictionBase):
16969a53589Sjeremylt    """Ceed Strided ElemRestriction: strided restriction from local vectors to elements."""
170c8e769f6Sjeremylt
171c8e769f6Sjeremylt    # Constructor
172d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, lsize, strides):
173c8e769f6Sjeremylt        # CeedVector object
174c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
175c8e769f6Sjeremylt
176c8e769f6Sjeremylt        # Reference to Ceed
177c8e769f6Sjeremylt        self._ceed = ceed
178c8e769f6Sjeremylt
17969a53589Sjeremylt        # Setup the numpy array for the libCEED call
18069a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
18169a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
18269a53589Sjeremylt                                   strides.__array_interface__['data'][0])
18369a53589Sjeremylt
184c8e769f6Sjeremylt        # libCEED call
185*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0],
186*477729cfSJeremy L Thompson                                                        nelem, elemsize, ncomp,
187*477729cfSJeremy L Thompson                                                        lsize, strides_pointer,
188*477729cfSJeremy L Thompson                                                        self._pointer)
189*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
190c8e769f6Sjeremylt
191c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1927a7b0fa3SJed Brown
1937a7b0fa3SJed Brown
194c8e769f6Sjeremyltclass BlockedElemRestriction(_ElemRestrictionBase):
195c8e769f6Sjeremylt    """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements."""
196c8e769f6Sjeremylt
197c8e769f6Sjeremylt    # Constructor
198d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize,
199d979a051Sjeremylt                 offsets, memtype=MEM_HOST, cmode=COPY_VALUES):
200c8e769f6Sjeremylt        # CeedVector object
201c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
202c8e769f6Sjeremylt
203c8e769f6Sjeremylt        # Reference to Ceed
204c8e769f6Sjeremylt        self._ceed = ceed
205c8e769f6Sjeremylt
206c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
207d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
208d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
209d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
210c8e769f6Sjeremylt
211c8e769f6Sjeremylt        # libCEED call
212*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], nelem,
213d979a051Sjeremylt                                                        elemsize, blksize, ncomp,
214d979a051Sjeremylt                                                        compstride, lsize, memtype, cmode,
215d979a051Sjeremylt                                                        offsets_pointer, self._pointer)
216*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
217c8e769f6Sjeremylt
218c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
219c8e769f6Sjeremylt    @property
220c8e769f6Sjeremylt    def T(self):
221c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
222c8e769f6Sjeremylt
223c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
224c8e769f6Sjeremylt
225c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
226c8e769f6Sjeremylt    @property
227c8e769f6Sjeremylt    def transpose(self):
228c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
229c8e769f6Sjeremylt
230c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
231c8e769f6Sjeremylt
232c8e769f6Sjeremylt    # Apply CeedElemRestriction to single block
233a8d32208Sjeremylt    def apply_block(self, block, u, v, tmode=NOTRANSPOSE,
234c8e769f6Sjeremylt                    request=REQUEST_IMMEDIATE):
235c8e769f6Sjeremylt        """Restrict a local vector to a block of an element vector or apply its transpose.
236c8e769f6Sjeremylt
237c8e769f6Sjeremylt           Args:
238c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
239c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
240c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
241c8e769f6Sjeremylt             u: input vector
242c8e769f6Sjeremylt             v: output vector
243c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
244c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
245c8e769f6Sjeremylt
246c8e769f6Sjeremylt        # libCEED call
247*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode,
248*477729cfSJeremy L Thompson                                                     u._pointer[0], v._pointer[0],
249*477729cfSJeremy L Thompson                                                     request)
250*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
251c8e769f6Sjeremylt
252c8e769f6Sjeremylt# ------------------------------------------------------------------------------
2537a7b0fa3SJed Brown
2547a7b0fa3SJed Brown
25569a53589Sjeremyltclass BlockedStridedElemRestriction(BlockedElemRestriction):
25669a53589Sjeremylt    """Ceed Blocked Strided ElemRestriction: strided restriction from local vectors to elements."""
25769a53589Sjeremylt
25869a53589Sjeremylt    # Constructor
259d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, lsize, strides):
26069a53589Sjeremylt        # CeedVector object
26169a53589Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
26269a53589Sjeremylt
26369a53589Sjeremylt        # Reference to Ceed
26469a53589Sjeremylt        self._ceed = ceed
26569a53589Sjeremylt
26669a53589Sjeremylt        # Setup the numpy array for the libCEED call
26769a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
26869a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
26969a53589Sjeremylt                                   strides.__array_interface__['data'][0])
27069a53589Sjeremylt
27169a53589Sjeremylt        # libCEED call
272*477729cfSJeremy L Thompson        err_code = lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem,
273d979a051Sjeremylt                                                               elemsize, blksize, ncomp,
274d979a051Sjeremylt                                                               lsize, strides_pointer,
27569a53589Sjeremylt                                                               self._pointer)
276*477729cfSJeremy L Thompson        self._ceed._check_error(err_code)
27769a53589Sjeremylt
27869a53589Sjeremylt# ------------------------------------------------------------------------------
2797a7b0fa3SJed Brown
2807a7b0fa3SJed Brown
281c8e769f6Sjeremyltclass TransposeElemRestriction():
282c8e769f6Sjeremylt    """Ceed ElemRestriction: transpose restriction from elements to local vectors."""
283c8e769f6Sjeremylt
284c8e769f6Sjeremylt    # Attributes
285c8e769f6Sjeremylt    _elemrestriction = None
286c8e769f6Sjeremylt
287c8e769f6Sjeremylt    # Constructor
288c8e769f6Sjeremylt    def __init__(self, elemrestriction):
289c8e769f6Sjeremylt
290c8e769f6Sjeremylt        # Reference elemrestriction
291c8e769f6Sjeremylt        self._elemrestriction = elemrestriction
292c8e769f6Sjeremylt
293c8e769f6Sjeremylt    # Representation
294c8e769f6Sjeremylt    def __repr__(self):
2957a7b0fa3SJed Brown        return "<Transpose CeedElemRestriction instance at " + \
2967a7b0fa3SJed Brown            hex(id(self)) + ">"
297c8e769f6Sjeremylt
298c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
2997a7b0fa3SJed Brown
300a8d32208Sjeremylt    def apply(self, u, v, request=REQUEST_IMMEDIATE):
301c8e769f6Sjeremylt        """Restrict an element vector to a local vector.
302c8e769f6Sjeremylt
303c8e769f6Sjeremylt           Args:
304c8e769f6Sjeremylt             u: input vector
305c8e769f6Sjeremylt             v: output vector
306c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
307c8e769f6Sjeremylt
308c8e769f6Sjeremylt        # libCEED call
309a8d32208Sjeremylt        self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE)
310c8e769f6Sjeremylt
311c8e769f6Sjeremylt# ------------------------------------------------------------------------------
3127a7b0fa3SJed Brown
3137a7b0fa3SJed Brown
314c8e769f6Sjeremyltclass TransposeBlockedElemRestriction(TransposeElemRestriction):
315c8e769f6Sjeremylt    """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements
316c8e769f6Sjeremylt         to local vectors."""
317c8e769f6Sjeremylt
318c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
319a8d32208Sjeremylt    def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE):
320c8e769f6Sjeremylt        """Restrict a block of an element vector to a local vector.
321c8e769f6Sjeremylt
322c8e769f6Sjeremylt           Args:
323c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
324c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
325c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
326c8e769f6Sjeremylt             u: input vector
327c8e769f6Sjeremylt             v: output vector
328c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
329c8e769f6Sjeremylt
330c8e769f6Sjeremylt        # libCEED call
331a8d32208Sjeremylt        self._elemrestriction.apply_block(block, u, v, request=request,
332c8e769f6Sjeremylt                                          tmode=TRANSPOSE)
333c8e769f6Sjeremylt
334c8e769f6Sjeremylt# ------------------------------------------------------------------------------
335