xref: /libCEED/python/ceed_elemrestriction.py (revision d979a0510f6353f9b8b50621433d53955a3f350a)
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*d979a051Sjeremyltfrom .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
37c8e769f6Sjeremylt        lib.CeedElemRestrictionDestroy(self._pointer)
38c8e769f6Sjeremylt
39c8e769f6Sjeremylt    # Representation
40c8e769f6Sjeremylt    def __repr__(self):
41c8e769f6Sjeremylt        return "<CeedElemRestriction instance at " + hex(id(self)) + ">"
42c8e769f6Sjeremylt
43c8e769f6Sjeremylt    # String conversion for print() to stdout
44c8e769f6Sjeremylt    def __str__(self):
45c8e769f6Sjeremylt        """View an ElemRestriction via print()."""
46c8e769f6Sjeremylt
47c8e769f6Sjeremylt        # libCEED call
48c8e769f6Sjeremylt        with tempfile.NamedTemporaryFile() as key_file:
49c8e769f6Sjeremylt            with open(key_file.name, 'r+') as stream_file:
50c8e769f6Sjeremylt                stream = ffi.cast("FILE *", stream_file)
51c8e769f6Sjeremylt
52c8e769f6Sjeremylt                lib.CeedElemRestrictionView(self._pointer[0], stream)
53c8e769f6Sjeremylt
54c8e769f6Sjeremylt                stream_file.seek(0)
55c8e769f6Sjeremylt                out_string = stream_file.read()
56c8e769f6Sjeremylt
57c8e769f6Sjeremylt        return out_string
58c8e769f6Sjeremylt
59c8e769f6Sjeremylt    # Apply CeedElemRestriction
60a8d32208Sjeremylt    def apply(self, u, v, tmode=NOTRANSPOSE, request=REQUEST_IMMEDIATE):
61c8e769f6Sjeremylt        """Restrict a local vector to an element vector or apply its transpose.
62c8e769f6Sjeremylt
63c8e769f6Sjeremylt           Args:
64c8e769f6Sjeremylt             u: input vector
65c8e769f6Sjeremylt             v: output vector
66c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
67c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
68c8e769f6Sjeremylt
69c8e769f6Sjeremylt        # libCEED call
70a8d32208Sjeremylt        lib.CeedElemRestrictionApply(self._pointer[0], tmode, u._pointer[0],
71a8d32208Sjeremylt                                     v._pointer[0], request)
72c8e769f6Sjeremylt
73c8e769f6Sjeremylt    # Transpose an ElemRestriction
74c8e769f6Sjeremylt    @property
75c8e769f6Sjeremylt    def T(self):
76c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
77c8e769f6Sjeremylt
78c8e769f6Sjeremylt        return TransposeElemRestriction(self)
79c8e769f6Sjeremylt
80c8e769f6Sjeremylt    # Transpose an ElemRestriction
81c8e769f6Sjeremylt    @property
82c8e769f6Sjeremylt    def transpose(self):
83c8e769f6Sjeremylt        """Transpose an ElemRestriction."""
84c8e769f6Sjeremylt
85c8e769f6Sjeremylt        return TransposeElemRestriction(self)
86c8e769f6Sjeremylt
87c8e769f6Sjeremylt    # Create restriction vectors
88c8e769f6Sjeremylt    def create_vector(self, createLvec=True, createEvec=True):
89c8e769f6Sjeremylt        """Create Vectors associated with an ElemRestriction.
90c8e769f6Sjeremylt
91c8e769f6Sjeremylt           Args:
92c8e769f6Sjeremylt             **createLvec: flag to create local vector, default True
93c8e769f6Sjeremylt             **createEvec: flag to create element vector, default True
94c8e769f6Sjeremylt
95c8e769f6Sjeremylt           Returns:
96c8e769f6Sjeremylt             [lvec, evec]: local vector and element vector, or None if flag set to false"""
97c8e769f6Sjeremylt
98c8e769f6Sjeremylt        # Vector pointers
99c8e769f6Sjeremylt        lvecPointer = ffi.new("CeedVector *") if createLvec else ffi.NULL
100c8e769f6Sjeremylt        evecPointer = ffi.new("CeedVector *") if createEvec else ffi.NULL
101c8e769f6Sjeremylt
102c8e769f6Sjeremylt        # libCEED call
103c8e769f6Sjeremylt        lib.CeedElemRestrictionCreateVector(self._pointer[0], lvecPointer,
104c8e769f6Sjeremylt                                            evecPointer)
105c8e769f6Sjeremylt
106c8e769f6Sjeremylt        # Return vectors
1077a7b0fa3SJed Brown        lvec = _VectorWrap(
1087a7b0fa3SJed Brown            self._ceed._pointer,
1097a7b0fa3SJed Brown            lvecPointer) if createLvec else None
1107a7b0fa3SJed Brown        evec = _VectorWrap(
1117a7b0fa3SJed Brown            self._ceed._pointer,
1127a7b0fa3SJed Brown            evecPointer) if createEvec else None
113c8e769f6Sjeremylt
114c8e769f6Sjeremylt        # Return
115c8e769f6Sjeremylt        return [lvec, evec]
116c8e769f6Sjeremylt
117c8e769f6Sjeremylt    # Get ElemRestriction multiplicity
118a8d32208Sjeremylt    def get_multiplicity(self):
119c8e769f6Sjeremylt        """Get the multiplicity of nodes in an ElemRestriction.
120c8e769f6Sjeremylt
121c8e769f6Sjeremylt           Returns:
122c8e769f6Sjeremylt             mult: local vector containing multiplicity of nodes in ElemRestriction"""
123c8e769f6Sjeremylt
124c8e769f6Sjeremylt        # Create mult vector
125c8e769f6Sjeremylt        [mult, evec] = self.create_vector(createEvec=False)
126c8e769f6Sjeremylt        mult.set_value(0)
127c8e769f6Sjeremylt
128c8e769f6Sjeremylt        # libCEED call
1297a7b0fa3SJed Brown        lib.CeedElemRestrictionGetMultiplicity(
1307a7b0fa3SJed Brown            self._pointer[0], mult._pointer[0])
131c8e769f6Sjeremylt
132c8e769f6Sjeremylt        # Return
133c8e769f6Sjeremylt        return mult
134c8e769f6Sjeremylt
135c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1367a7b0fa3SJed Brown
1377a7b0fa3SJed Brown
138c8e769f6Sjeremyltclass ElemRestriction(_ElemRestrictionBase):
139c8e769f6Sjeremylt    """Ceed ElemRestriction: restriction from local vectors to elements."""
140c8e769f6Sjeremylt
141c8e769f6Sjeremylt    # Constructor
142*d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets,
143*d979a051Sjeremylt                 memtype=MEM_HOST, cmode=COPY_VALUES):
144c8e769f6Sjeremylt        # CeedVector object
145c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
146c8e769f6Sjeremylt
147c8e769f6Sjeremylt        # Reference to Ceed
148c8e769f6Sjeremylt        self._ceed = ceed
149c8e769f6Sjeremylt
150c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
151*d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
152*d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
153*d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
154c8e769f6Sjeremylt
155c8e769f6Sjeremylt        # libCEED call
156*d979a051Sjeremylt        lib.CeedElemRestrictionCreate(self._ceed._pointer[0], nelem, elemsize,
157*d979a051Sjeremylt                                      ncomp, compstride, lsize, memtype, cmode,
158*d979a051Sjeremylt                                      offsets_pointer, self._pointer)
159c8e769f6Sjeremylt
160c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1617a7b0fa3SJed Brown
1627a7b0fa3SJed Brown
16369a53589Sjeremyltclass StridedElemRestriction(_ElemRestrictionBase):
16469a53589Sjeremylt    """Ceed Strided ElemRestriction: strided restriction from local vectors to elements."""
165c8e769f6Sjeremylt
166c8e769f6Sjeremylt    # Constructor
167*d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, ncomp, lsize, strides):
168c8e769f6Sjeremylt        # CeedVector object
169c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
170c8e769f6Sjeremylt
171c8e769f6Sjeremylt        # Reference to Ceed
172c8e769f6Sjeremylt        self._ceed = ceed
173c8e769f6Sjeremylt
17469a53589Sjeremylt        # Setup the numpy array for the libCEED call
17569a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
17669a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
17769a53589Sjeremylt                                   strides.__array_interface__['data'][0])
17869a53589Sjeremylt
179c8e769f6Sjeremylt        # libCEED call
18069a53589Sjeremylt        lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0], nelem,
181*d979a051Sjeremylt                                             elemsize, ncomp, lsize,
18269a53589Sjeremylt                                             strides_pointer, self._pointer)
183c8e769f6Sjeremylt
184c8e769f6Sjeremylt# ------------------------------------------------------------------------------
1857a7b0fa3SJed Brown
1867a7b0fa3SJed Brown
187c8e769f6Sjeremyltclass BlockedElemRestriction(_ElemRestrictionBase):
188c8e769f6Sjeremylt    """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements."""
189c8e769f6Sjeremylt
190c8e769f6Sjeremylt    # Constructor
191*d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize,
192*d979a051Sjeremylt                 offsets, memtype=MEM_HOST, cmode=COPY_VALUES):
193c8e769f6Sjeremylt        # CeedVector object
194c8e769f6Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
195c8e769f6Sjeremylt
196c8e769f6Sjeremylt        # Reference to Ceed
197c8e769f6Sjeremylt        self._ceed = ceed
198c8e769f6Sjeremylt
199c8e769f6Sjeremylt        # Setup the numpy array for the libCEED call
200*d979a051Sjeremylt        offsets_pointer = ffi.new("const CeedInt *")
201*d979a051Sjeremylt        offsets_pointer = ffi.cast("const CeedInt *",
202*d979a051Sjeremylt                                   offsets.__array_interface__['data'][0])
203c8e769f6Sjeremylt
204c8e769f6Sjeremylt        # libCEED call
205*d979a051Sjeremylt        lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], nelem,
206*d979a051Sjeremylt                                             elemsize, blksize, ncomp,
207*d979a051Sjeremylt                                             compstride, lsize, memtype, cmode,
208*d979a051Sjeremylt                                             offsets_pointer, self._pointer)
209c8e769f6Sjeremylt
210c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
211c8e769f6Sjeremylt    @property
212c8e769f6Sjeremylt    def T(self):
213c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
214c8e769f6Sjeremylt
215c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
216c8e769f6Sjeremylt
217c8e769f6Sjeremylt    # Transpose a Blocked ElemRestriction
218c8e769f6Sjeremylt    @property
219c8e769f6Sjeremylt    def transpose(self):
220c8e769f6Sjeremylt        """Transpose a BlockedElemRestriction."""
221c8e769f6Sjeremylt
222c8e769f6Sjeremylt        return TransposeBlockedElemRestriction(self)
223c8e769f6Sjeremylt
224c8e769f6Sjeremylt    # Apply CeedElemRestriction to single block
225a8d32208Sjeremylt    def apply_block(self, block, u, v, tmode=NOTRANSPOSE,
226c8e769f6Sjeremylt                    request=REQUEST_IMMEDIATE):
227c8e769f6Sjeremylt        """Restrict a local vector to a block of an element vector or apply its transpose.
228c8e769f6Sjeremylt
229c8e769f6Sjeremylt           Args:
230c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
231c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
232c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
233c8e769f6Sjeremylt             u: input vector
234c8e769f6Sjeremylt             v: output vector
235c8e769f6Sjeremylt             **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
236c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
237c8e769f6Sjeremylt
238c8e769f6Sjeremylt        # libCEED call
239c8e769f6Sjeremylt        lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode,
240a8d32208Sjeremylt                                          u._pointer[0], v._pointer[0], request)
241c8e769f6Sjeremylt
242c8e769f6Sjeremylt# ------------------------------------------------------------------------------
2437a7b0fa3SJed Brown
2447a7b0fa3SJed Brown
24569a53589Sjeremyltclass BlockedStridedElemRestriction(BlockedElemRestriction):
24669a53589Sjeremylt    """Ceed Blocked Strided ElemRestriction: strided restriction from local vectors to elements."""
24769a53589Sjeremylt
24869a53589Sjeremylt    # Constructor
249*d979a051Sjeremylt    def __init__(self, ceed, nelem, elemsize, blksize, ncomp, lsize, strides):
25069a53589Sjeremylt        # CeedVector object
25169a53589Sjeremylt        self._pointer = ffi.new("CeedElemRestriction *")
25269a53589Sjeremylt
25369a53589Sjeremylt        # Reference to Ceed
25469a53589Sjeremylt        self._ceed = ceed
25569a53589Sjeremylt
25669a53589Sjeremylt        # Setup the numpy array for the libCEED call
25769a53589Sjeremylt        strides_pointer = ffi.new("const CeedInt *")
25869a53589Sjeremylt        strides_pointer = ffi.cast("const CeedInt *",
25969a53589Sjeremylt                                   strides.__array_interface__['data'][0])
26069a53589Sjeremylt
26169a53589Sjeremylt        # libCEED call
26269a53589Sjeremylt        lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem,
263*d979a051Sjeremylt                                                    elemsize, blksize, ncomp,
264*d979a051Sjeremylt                                                    lsize, strides_pointer,
26569a53589Sjeremylt                                                    self._pointer)
26669a53589Sjeremylt
26769a53589Sjeremylt# ------------------------------------------------------------------------------
2687a7b0fa3SJed Brown
2697a7b0fa3SJed Brown
270c8e769f6Sjeremyltclass TransposeElemRestriction():
271c8e769f6Sjeremylt    """Ceed ElemRestriction: transpose restriction from elements to local vectors."""
272c8e769f6Sjeremylt
273c8e769f6Sjeremylt    # Attributes
274c8e769f6Sjeremylt    _elemrestriction = None
275c8e769f6Sjeremylt
276c8e769f6Sjeremylt    # Constructor
277c8e769f6Sjeremylt    def __init__(self, elemrestriction):
278c8e769f6Sjeremylt
279c8e769f6Sjeremylt        # Reference elemrestriction
280c8e769f6Sjeremylt        self._elemrestriction = elemrestriction
281c8e769f6Sjeremylt
282c8e769f6Sjeremylt    # Representation
283c8e769f6Sjeremylt    def __repr__(self):
2847a7b0fa3SJed Brown        return "<Transpose CeedElemRestriction instance at " + \
2857a7b0fa3SJed Brown            hex(id(self)) + ">"
286c8e769f6Sjeremylt
287c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
2887a7b0fa3SJed Brown
289a8d32208Sjeremylt    def apply(self, u, v, request=REQUEST_IMMEDIATE):
290c8e769f6Sjeremylt        """Restrict an element vector to a local vector.
291c8e769f6Sjeremylt
292c8e769f6Sjeremylt           Args:
293c8e769f6Sjeremylt             u: input vector
294c8e769f6Sjeremylt             v: output vector
295c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
296c8e769f6Sjeremylt
297c8e769f6Sjeremylt        # libCEED call
298a8d32208Sjeremylt        self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE)
299c8e769f6Sjeremylt
300c8e769f6Sjeremylt# ------------------------------------------------------------------------------
3017a7b0fa3SJed Brown
3027a7b0fa3SJed Brown
303c8e769f6Sjeremyltclass TransposeBlockedElemRestriction(TransposeElemRestriction):
304c8e769f6Sjeremylt    """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements
305c8e769f6Sjeremylt         to local vectors."""
306c8e769f6Sjeremylt
307c8e769f6Sjeremylt    # Apply Transpose CeedElemRestriction
308a8d32208Sjeremylt    def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE):
309c8e769f6Sjeremylt        """Restrict a block of an element vector to a local vector.
310c8e769f6Sjeremylt
311c8e769f6Sjeremylt           Args:
312c8e769f6Sjeremylt             block: block number to restrict to/from, i.e. block=0 will handle
313c8e769f6Sjeremylt                      elements [0 : blksize] and block=3 will handle elements
314c8e769f6Sjeremylt                      [3*blksize : 4*blksize]
315c8e769f6Sjeremylt             u: input vector
316c8e769f6Sjeremylt             v: output vector
317c8e769f6Sjeremylt             **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
318c8e769f6Sjeremylt
319c8e769f6Sjeremylt        # libCEED call
320a8d32208Sjeremylt        self._elemrestriction.apply_block(block, u, v, request=request,
321c8e769f6Sjeremylt                                          tmode=TRANSPOSE)
322c8e769f6Sjeremylt
323c8e769f6Sjeremylt# ------------------------------------------------------------------------------
324