xref: /libCEED/python/ceed_elemrestriction.py (revision 69a53589b8c3c82ae8bb5be56f1100c7d3468573)
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
2161dbc9d2Sjeremyltfrom .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, MEM_HOST, COPY_VALUES, TRANSPOSE, NOTRANSPOSE, INTERLACED, NONINTERLACED
22c8e769f6Sjeremyltfrom .ceed_vector import _VectorWrap
23c8e769f6Sjeremylt
24c8e769f6Sjeremylt# ------------------------------------------------------------------------------
25c8e769f6Sjeremyltclass _ElemRestrictionBase(ABC):
26c8e769f6Sjeremylt  """Ceed ElemRestriction: restriction from local vectors to elements."""
27c8e769f6Sjeremylt
28c8e769f6Sjeremylt  # Attributes
29c8e769f6Sjeremylt  _ceed = ffi.NULL
30c8e769f6Sjeremylt  _pointer = ffi.NULL
31c8e769f6Sjeremylt
32c8e769f6Sjeremylt  # Destructor
33c8e769f6Sjeremylt  def __del__(self):
34c8e769f6Sjeremylt    # libCEED call
35c8e769f6Sjeremylt    lib.CeedElemRestrictionDestroy(self._pointer)
36c8e769f6Sjeremylt
37c8e769f6Sjeremylt  # Representation
38c8e769f6Sjeremylt  def __repr__(self):
39c8e769f6Sjeremylt    return "<CeedElemRestriction instance at " + hex(id(self)) + ">"
40c8e769f6Sjeremylt
41c8e769f6Sjeremylt  # String conversion for print() to stdout
42c8e769f6Sjeremylt  def __str__(self):
43c8e769f6Sjeremylt    """View an ElemRestriction via print()."""
44c8e769f6Sjeremylt
45c8e769f6Sjeremylt    # libCEED call
46c8e769f6Sjeremylt    with tempfile.NamedTemporaryFile() as key_file:
47c8e769f6Sjeremylt      with open(key_file.name, 'r+') as stream_file:
48c8e769f6Sjeremylt        stream = ffi.cast("FILE *", stream_file)
49c8e769f6Sjeremylt
50c8e769f6Sjeremylt        lib.CeedElemRestrictionView(self._pointer[0], stream)
51c8e769f6Sjeremylt
52c8e769f6Sjeremylt        stream_file.seek(0)
53c8e769f6Sjeremylt        out_string = stream_file.read()
54c8e769f6Sjeremylt
55c8e769f6Sjeremylt    return out_string
56c8e769f6Sjeremylt
57c8e769f6Sjeremylt  # Apply CeedElemRestriction
58a8d32208Sjeremylt  def apply(self, u, v, tmode=NOTRANSPOSE, request=REQUEST_IMMEDIATE):
59c8e769f6Sjeremylt    """Restrict a local vector to an element vector or apply its transpose.
60c8e769f6Sjeremylt
61c8e769f6Sjeremylt       Args:
62c8e769f6Sjeremylt         u: input vector
63c8e769f6Sjeremylt         v: output vector
64c8e769f6Sjeremylt         **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
65c8e769f6Sjeremylt         **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
66c8e769f6Sjeremylt
67c8e769f6Sjeremylt    # libCEED call
68a8d32208Sjeremylt    lib.CeedElemRestrictionApply(self._pointer[0], tmode, u._pointer[0],
69a8d32208Sjeremylt                                 v._pointer[0], request)
70c8e769f6Sjeremylt
71c8e769f6Sjeremylt  # Transpose an ElemRestriction
72c8e769f6Sjeremylt  @property
73c8e769f6Sjeremylt  def T(self):
74c8e769f6Sjeremylt    """Transpose an ElemRestriction."""
75c8e769f6Sjeremylt
76c8e769f6Sjeremylt    return TransposeElemRestriction(self)
77c8e769f6Sjeremylt
78c8e769f6Sjeremylt  # Transpose an ElemRestriction
79c8e769f6Sjeremylt  @property
80c8e769f6Sjeremylt  def transpose(self):
81c8e769f6Sjeremylt    """Transpose an ElemRestriction."""
82c8e769f6Sjeremylt
83c8e769f6Sjeremylt    return TransposeElemRestriction(self)
84c8e769f6Sjeremylt
85c8e769f6Sjeremylt  # Create restriction vectors
86c8e769f6Sjeremylt  def create_vector(self, createLvec = True, createEvec = True):
87c8e769f6Sjeremylt    """Create Vectors associated with an ElemRestriction.
88c8e769f6Sjeremylt
89c8e769f6Sjeremylt       Args:
90c8e769f6Sjeremylt         **createLvec: flag to create local vector, default True
91c8e769f6Sjeremylt         **createEvec: flag to create element vector, default True
92c8e769f6Sjeremylt
93c8e769f6Sjeremylt       Returns:
94c8e769f6Sjeremylt         [lvec, evec]: local vector and element vector, or None if flag set to false"""
95c8e769f6Sjeremylt
96c8e769f6Sjeremylt    # Vector pointers
97c8e769f6Sjeremylt    lvecPointer = ffi.new("CeedVector *") if createLvec else ffi.NULL
98c8e769f6Sjeremylt    evecPointer = ffi.new("CeedVector *") if createEvec else ffi.NULL
99c8e769f6Sjeremylt
100c8e769f6Sjeremylt    # libCEED call
101c8e769f6Sjeremylt    lib.CeedElemRestrictionCreateVector(self._pointer[0], lvecPointer,
102c8e769f6Sjeremylt                                        evecPointer)
103c8e769f6Sjeremylt
104c8e769f6Sjeremylt    # Return vectors
105c8e769f6Sjeremylt    lvec = _VectorWrap(self._ceed._pointer, lvecPointer) if createLvec else None
106c8e769f6Sjeremylt    evec = _VectorWrap(self._ceed._pointer, evecPointer) if createEvec else None
107c8e769f6Sjeremylt
108c8e769f6Sjeremylt    # Return
109c8e769f6Sjeremylt    return [lvec, evec]
110c8e769f6Sjeremylt
111c8e769f6Sjeremylt  # Get ElemRestriction multiplicity
112a8d32208Sjeremylt  def get_multiplicity(self):
113c8e769f6Sjeremylt    """Get the multiplicity of nodes in an ElemRestriction.
114c8e769f6Sjeremylt
115c8e769f6Sjeremylt       Returns:
116c8e769f6Sjeremylt         mult: local vector containing multiplicity of nodes in ElemRestriction"""
117c8e769f6Sjeremylt
118c8e769f6Sjeremylt    # Create mult vector
119c8e769f6Sjeremylt    [mult, evec] = self.create_vector(createEvec = False)
120c8e769f6Sjeremylt    mult.set_value(0)
121c8e769f6Sjeremylt
122c8e769f6Sjeremylt    # libCEED call
123a8d32208Sjeremylt    lib.CeedElemRestrictionGetMultiplicity(self._pointer[0], mult._pointer[0])
124c8e769f6Sjeremylt
125c8e769f6Sjeremylt    # Return
126c8e769f6Sjeremylt    return mult
127c8e769f6Sjeremylt
128c8e769f6Sjeremylt# ------------------------------------------------------------------------------
129c8e769f6Sjeremyltclass ElemRestriction(_ElemRestrictionBase):
130c8e769f6Sjeremylt  """Ceed ElemRestriction: restriction from local vectors to elements."""
131c8e769f6Sjeremylt
132c8e769f6Sjeremylt  # Constructor
133c8e769f6Sjeremylt  def __init__(self, ceed, nelem, elemsize, nnodes, ncomp, indices,
13461dbc9d2Sjeremylt               memtype=MEM_HOST, cmode=COPY_VALUES, imode=NONINTERLACED):
135c8e769f6Sjeremylt    # CeedVector object
136c8e769f6Sjeremylt    self._pointer = ffi.new("CeedElemRestriction *")
137c8e769f6Sjeremylt
138c8e769f6Sjeremylt    # Reference to Ceed
139c8e769f6Sjeremylt    self._ceed = ceed
140c8e769f6Sjeremylt
141c8e769f6Sjeremylt    # Setup the numpy array for the libCEED call
142c8e769f6Sjeremylt    indices_pointer = ffi.new("const CeedInt *")
143c8e769f6Sjeremylt    indices_pointer = ffi.cast("const CeedInt *",
144c8e769f6Sjeremylt                               indices.__array_interface__['data'][0])
145c8e769f6Sjeremylt
146c8e769f6Sjeremylt    # libCEED call
14761dbc9d2Sjeremylt    lib.CeedElemRestrictionCreate(self._ceed._pointer[0], imode, nelem,
148a8d32208Sjeremylt                                  elemsize, nnodes, ncomp, memtype, cmode,
149c8e769f6Sjeremylt                                  indices_pointer, self._pointer)
150c8e769f6Sjeremylt
151c8e769f6Sjeremylt# ------------------------------------------------------------------------------
152*69a53589Sjeremyltclass StridedElemRestriction(_ElemRestrictionBase):
153*69a53589Sjeremylt  """Ceed Strided ElemRestriction: strided restriction from local vectors to elements."""
154c8e769f6Sjeremylt
155c8e769f6Sjeremylt  # Constructor
156*69a53589Sjeremylt  def __init__(self, ceed, nelem, elemsize, nnodes, ncomp, strides):
157c8e769f6Sjeremylt    # CeedVector object
158c8e769f6Sjeremylt    self._pointer = ffi.new("CeedElemRestriction *")
159c8e769f6Sjeremylt
160c8e769f6Sjeremylt    # Reference to Ceed
161c8e769f6Sjeremylt    self._ceed = ceed
162c8e769f6Sjeremylt
163*69a53589Sjeremylt    # Setup the numpy array for the libCEED call
164*69a53589Sjeremylt    strides_pointer = ffi.new("const CeedInt *")
165*69a53589Sjeremylt    strides_pointer = ffi.cast("const CeedInt *",
166*69a53589Sjeremylt                               strides.__array_interface__['data'][0])
167*69a53589Sjeremylt
168c8e769f6Sjeremylt    # libCEED call
169*69a53589Sjeremylt    lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0], nelem,
170c8e769f6Sjeremylt                                         elemsize, nnodes, ncomp,
171*69a53589Sjeremylt                                         strides_pointer, self._pointer)
172c8e769f6Sjeremylt
173c8e769f6Sjeremylt# ------------------------------------------------------------------------------
174c8e769f6Sjeremyltclass BlockedElemRestriction(_ElemRestrictionBase):
175c8e769f6Sjeremylt  """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements."""
176c8e769f6Sjeremylt
177c8e769f6Sjeremylt  # Constructor
178c8e769f6Sjeremylt  def __init__(self, ceed, nelem, elemsize, blksize, nnodes, ncomp, indices,
17961dbc9d2Sjeremylt               memtype=MEM_HOST, cmode=COPY_VALUES, imode=NONINTERLACED):
180c8e769f6Sjeremylt    # CeedVector object
181c8e769f6Sjeremylt    self._pointer = ffi.new("CeedElemRestriction *")
182c8e769f6Sjeremylt
183c8e769f6Sjeremylt    # Reference to Ceed
184c8e769f6Sjeremylt    self._ceed = ceed
185c8e769f6Sjeremylt
186c8e769f6Sjeremylt    # Setup the numpy array for the libCEED call
187c8e769f6Sjeremylt    indices_pointer = ffi.new("const CeedInt *")
188c8e769f6Sjeremylt    indices_pointer = ffi.cast("const CeedInt *",
189c8e769f6Sjeremylt                               indices.__array_interface__['data'][0])
190c8e769f6Sjeremylt
191c8e769f6Sjeremylt    # libCEED call
19261dbc9d2Sjeremylt    lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], imode, nelem,
193c8e769f6Sjeremylt                                         elemsize, blksize, nnodes, ncomp,
194c8e769f6Sjeremylt                                         memtype, cmode, indices_pointer,
195c8e769f6Sjeremylt                                         self._pointer)
196c8e769f6Sjeremylt
197c8e769f6Sjeremylt  # Transpose a Blocked ElemRestriction
198c8e769f6Sjeremylt  @property
199c8e769f6Sjeremylt  def T(self):
200c8e769f6Sjeremylt    """Transpose a BlockedElemRestriction."""
201c8e769f6Sjeremylt
202c8e769f6Sjeremylt    return TransposeBlockedElemRestriction(self)
203c8e769f6Sjeremylt
204c8e769f6Sjeremylt  # Transpose a Blocked ElemRestriction
205c8e769f6Sjeremylt  @property
206c8e769f6Sjeremylt  def transpose(self):
207c8e769f6Sjeremylt    """Transpose a BlockedElemRestriction."""
208c8e769f6Sjeremylt
209c8e769f6Sjeremylt    return TransposeBlockedElemRestriction(self)
210c8e769f6Sjeremylt
211c8e769f6Sjeremylt  # Apply CeedElemRestriction to single block
212a8d32208Sjeremylt  def apply_block(self, block, u, v, tmode=NOTRANSPOSE,
213c8e769f6Sjeremylt                  request=REQUEST_IMMEDIATE):
214c8e769f6Sjeremylt    """Restrict a local vector to a block of an element vector or apply its transpose.
215c8e769f6Sjeremylt
216c8e769f6Sjeremylt       Args:
217c8e769f6Sjeremylt         block: block number to restrict to/from, i.e. block=0 will handle
218c8e769f6Sjeremylt                  elements [0 : blksize] and block=3 will handle elements
219c8e769f6Sjeremylt                  [3*blksize : 4*blksize]
220c8e769f6Sjeremylt         u: input vector
221c8e769f6Sjeremylt         v: output vector
222c8e769f6Sjeremylt         **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE
223c8e769f6Sjeremylt         **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
224c8e769f6Sjeremylt
225c8e769f6Sjeremylt    # libCEED call
226c8e769f6Sjeremylt    lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode,
227a8d32208Sjeremylt                                      u._pointer[0], v._pointer[0], request)
228c8e769f6Sjeremylt
229c8e769f6Sjeremylt# ------------------------------------------------------------------------------
230*69a53589Sjeremyltclass BlockedStridedElemRestriction(BlockedElemRestriction):
231*69a53589Sjeremylt  """Ceed Blocked Strided ElemRestriction: strided restriction from local vectors to elements."""
232*69a53589Sjeremylt
233*69a53589Sjeremylt  # Constructor
234*69a53589Sjeremylt  def __init__(self, ceed, nelem, elemsize, blksize, nnodes, ncomp, strides):
235*69a53589Sjeremylt    # CeedVector object
236*69a53589Sjeremylt    self._pointer = ffi.new("CeedElemRestriction *")
237*69a53589Sjeremylt
238*69a53589Sjeremylt    # Reference to Ceed
239*69a53589Sjeremylt    self._ceed = ceed
240*69a53589Sjeremylt
241*69a53589Sjeremylt    # Setup the numpy array for the libCEED call
242*69a53589Sjeremylt    strides_pointer = ffi.new("const CeedInt *")
243*69a53589Sjeremylt    strides_pointer = ffi.cast("const CeedInt *",
244*69a53589Sjeremylt                               strides.__array_interface__['data'][0])
245*69a53589Sjeremylt
246*69a53589Sjeremylt    # libCEED call
247*69a53589Sjeremylt    lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem,
248*69a53589Sjeremylt                                                elemsize, blksize, nnodes,
249*69a53589Sjeremylt                                                ncomp, strides_pointer,
250*69a53589Sjeremylt                                                self._pointer)
251*69a53589Sjeremylt
252*69a53589Sjeremylt# ------------------------------------------------------------------------------
253c8e769f6Sjeremyltclass TransposeElemRestriction():
254c8e769f6Sjeremylt  """Ceed ElemRestriction: transpose restriction from elements to local vectors."""
255c8e769f6Sjeremylt
256c8e769f6Sjeremylt  # Attributes
257c8e769f6Sjeremylt  _elemrestriction = None
258c8e769f6Sjeremylt
259c8e769f6Sjeremylt  # Constructor
260c8e769f6Sjeremylt  def __init__(self, elemrestriction):
261c8e769f6Sjeremylt
262c8e769f6Sjeremylt    # Reference elemrestriction
263c8e769f6Sjeremylt    self._elemrestriction = elemrestriction
264c8e769f6Sjeremylt
265c8e769f6Sjeremylt  # Representation
266c8e769f6Sjeremylt  def __repr__(self):
267c8e769f6Sjeremylt    return "<Transpose CeedElemRestriction instance at " + hex(id(self)) + ">"
268c8e769f6Sjeremylt
269c8e769f6Sjeremylt
270c8e769f6Sjeremylt  # Apply Transpose CeedElemRestriction
271a8d32208Sjeremylt  def apply(self, u, v, request=REQUEST_IMMEDIATE):
272c8e769f6Sjeremylt    """Restrict an element vector to a local vector.
273c8e769f6Sjeremylt
274c8e769f6Sjeremylt       Args:
275c8e769f6Sjeremylt         u: input vector
276c8e769f6Sjeremylt         v: output vector
277c8e769f6Sjeremylt         **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
278c8e769f6Sjeremylt
279c8e769f6Sjeremylt    # libCEED call
280a8d32208Sjeremylt    self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE)
281c8e769f6Sjeremylt
282c8e769f6Sjeremylt# ------------------------------------------------------------------------------
283c8e769f6Sjeremyltclass TransposeBlockedElemRestriction(TransposeElemRestriction):
284c8e769f6Sjeremylt  """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements
285c8e769f6Sjeremylt       to local vectors."""
286c8e769f6Sjeremylt
287c8e769f6Sjeremylt  # Apply Transpose CeedElemRestriction
288a8d32208Sjeremylt  def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE):
289c8e769f6Sjeremylt    """Restrict a block of an element vector to a local vector.
290c8e769f6Sjeremylt
291c8e769f6Sjeremylt       Args:
292c8e769f6Sjeremylt         block: block number to restrict to/from, i.e. block=0 will handle
293c8e769f6Sjeremylt                  elements [0 : blksize] and block=3 will handle elements
294c8e769f6Sjeremylt                  [3*blksize : 4*blksize]
295c8e769f6Sjeremylt         u: input vector
296c8e769f6Sjeremylt         v: output vector
297c8e769f6Sjeremylt         **request: Ceed request, default CEED_REQUEST_IMMEDIATE"""
298c8e769f6Sjeremylt
299c8e769f6Sjeremylt    # libCEED call
300a8d32208Sjeremylt    self._elemrestriction.apply_block(block, u, v, request=request,
301c8e769f6Sjeremylt                                      tmode=TRANSPOSE)
302c8e769f6Sjeremylt
303c8e769f6Sjeremylt# ------------------------------------------------------------------------------
304