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*61dbc9d2Sjeremyltfrom .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, 134*61dbc9d2Sjeremylt 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 147*61dbc9d2Sjeremylt lib.CeedElemRestrictionCreate(self._ceed._pointer[0], imode, nelem, 148a8d32208Sjeremylt elemsize, nnodes, ncomp, memtype, cmode, 149c8e769f6Sjeremylt indices_pointer, self._pointer) 150c8e769f6Sjeremylt 151c8e769f6Sjeremylt# ------------------------------------------------------------------------------ 152c8e769f6Sjeremyltclass IdentityElemRestriction(_ElemRestrictionBase): 153c8e769f6Sjeremylt """Ceed Identity ElemRestriction: identity restriction from local vectors to elements.""" 154c8e769f6Sjeremylt 155c8e769f6Sjeremylt # Constructor 156*61dbc9d2Sjeremylt def __init__(self, ceed, nelem, elemsize, nnodes, ncomp, imode=NONINTERLACED): 157c8e769f6Sjeremylt # CeedVector object 158c8e769f6Sjeremylt self._pointer = ffi.new("CeedElemRestriction *") 159c8e769f6Sjeremylt 160c8e769f6Sjeremylt # Reference to Ceed 161c8e769f6Sjeremylt self._ceed = ceed 162c8e769f6Sjeremylt 163c8e769f6Sjeremylt # libCEED call 164*61dbc9d2Sjeremylt lib.CeedElemRestrictionCreateIdentity(self._ceed._pointer[0], imode, nelem, 165c8e769f6Sjeremylt elemsize, nnodes, ncomp, 166c8e769f6Sjeremylt self._pointer) 167c8e769f6Sjeremylt 168c8e769f6Sjeremylt# ------------------------------------------------------------------------------ 169c8e769f6Sjeremyltclass BlockedElemRestriction(_ElemRestrictionBase): 170c8e769f6Sjeremylt """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements.""" 171c8e769f6Sjeremylt 172c8e769f6Sjeremylt # Constructor 173c8e769f6Sjeremylt def __init__(self, ceed, nelem, elemsize, blksize, nnodes, ncomp, indices, 174*61dbc9d2Sjeremylt memtype=MEM_HOST, cmode=COPY_VALUES, imode=NONINTERLACED): 175c8e769f6Sjeremylt # CeedVector object 176c8e769f6Sjeremylt self._pointer = ffi.new("CeedElemRestriction *") 177c8e769f6Sjeremylt 178c8e769f6Sjeremylt # Reference to Ceed 179c8e769f6Sjeremylt self._ceed = ceed 180c8e769f6Sjeremylt 181c8e769f6Sjeremylt # Setup the numpy array for the libCEED call 182c8e769f6Sjeremylt indices_pointer = ffi.new("const CeedInt *") 183c8e769f6Sjeremylt indices_pointer = ffi.cast("const CeedInt *", 184c8e769f6Sjeremylt indices.__array_interface__['data'][0]) 185c8e769f6Sjeremylt 186c8e769f6Sjeremylt # libCEED call 187*61dbc9d2Sjeremylt lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], imode, nelem, 188c8e769f6Sjeremylt elemsize, blksize, nnodes, ncomp, 189c8e769f6Sjeremylt memtype, cmode, indices_pointer, 190c8e769f6Sjeremylt self._pointer) 191c8e769f6Sjeremylt 192c8e769f6Sjeremylt # Transpose a Blocked ElemRestriction 193c8e769f6Sjeremylt @property 194c8e769f6Sjeremylt def T(self): 195c8e769f6Sjeremylt """Transpose a BlockedElemRestriction.""" 196c8e769f6Sjeremylt 197c8e769f6Sjeremylt return TransposeBlockedElemRestriction(self) 198c8e769f6Sjeremylt 199c8e769f6Sjeremylt # Transpose a Blocked ElemRestriction 200c8e769f6Sjeremylt @property 201c8e769f6Sjeremylt def transpose(self): 202c8e769f6Sjeremylt """Transpose a BlockedElemRestriction.""" 203c8e769f6Sjeremylt 204c8e769f6Sjeremylt return TransposeBlockedElemRestriction(self) 205c8e769f6Sjeremylt 206c8e769f6Sjeremylt # Apply CeedElemRestriction to single block 207a8d32208Sjeremylt def apply_block(self, block, u, v, tmode=NOTRANSPOSE, 208c8e769f6Sjeremylt request=REQUEST_IMMEDIATE): 209c8e769f6Sjeremylt """Restrict a local vector to a block of an element vector or apply its transpose. 210c8e769f6Sjeremylt 211c8e769f6Sjeremylt Args: 212c8e769f6Sjeremylt block: block number to restrict to/from, i.e. block=0 will handle 213c8e769f6Sjeremylt elements [0 : blksize] and block=3 will handle elements 214c8e769f6Sjeremylt [3*blksize : 4*blksize] 215c8e769f6Sjeremylt u: input vector 216c8e769f6Sjeremylt v: output vector 217c8e769f6Sjeremylt **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE 218c8e769f6Sjeremylt **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 219c8e769f6Sjeremylt 220c8e769f6Sjeremylt # libCEED call 221c8e769f6Sjeremylt lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode, 222a8d32208Sjeremylt u._pointer[0], v._pointer[0], request) 223c8e769f6Sjeremylt 224c8e769f6Sjeremylt# ------------------------------------------------------------------------------ 225c8e769f6Sjeremyltclass TransposeElemRestriction(): 226c8e769f6Sjeremylt """Ceed ElemRestriction: transpose restriction from elements to local vectors.""" 227c8e769f6Sjeremylt 228c8e769f6Sjeremylt # Attributes 229c8e769f6Sjeremylt _elemrestriction = None 230c8e769f6Sjeremylt 231c8e769f6Sjeremylt # Constructor 232c8e769f6Sjeremylt def __init__(self, elemrestriction): 233c8e769f6Sjeremylt 234c8e769f6Sjeremylt # Reference elemrestriction 235c8e769f6Sjeremylt self._elemrestriction = elemrestriction 236c8e769f6Sjeremylt 237c8e769f6Sjeremylt # Representation 238c8e769f6Sjeremylt def __repr__(self): 239c8e769f6Sjeremylt return "<Transpose CeedElemRestriction instance at " + hex(id(self)) + ">" 240c8e769f6Sjeremylt 241c8e769f6Sjeremylt 242c8e769f6Sjeremylt # Apply Transpose CeedElemRestriction 243a8d32208Sjeremylt def apply(self, u, v, request=REQUEST_IMMEDIATE): 244c8e769f6Sjeremylt """Restrict an element vector to a local vector. 245c8e769f6Sjeremylt 246c8e769f6Sjeremylt Args: 247c8e769f6Sjeremylt u: input vector 248c8e769f6Sjeremylt v: output vector 249c8e769f6Sjeremylt **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 250c8e769f6Sjeremylt 251c8e769f6Sjeremylt # libCEED call 252a8d32208Sjeremylt self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE) 253c8e769f6Sjeremylt 254c8e769f6Sjeremylt# ------------------------------------------------------------------------------ 255c8e769f6Sjeremyltclass TransposeBlockedElemRestriction(TransposeElemRestriction): 256c8e769f6Sjeremylt """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements 257c8e769f6Sjeremylt to local vectors.""" 258c8e769f6Sjeremylt 259c8e769f6Sjeremylt # Apply Transpose CeedElemRestriction 260a8d32208Sjeremylt def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE): 261c8e769f6Sjeremylt """Restrict a block of an element vector to a local vector. 262c8e769f6Sjeremylt 263c8e769f6Sjeremylt Args: 264c8e769f6Sjeremylt block: block number to restrict to/from, i.e. block=0 will handle 265c8e769f6Sjeremylt elements [0 : blksize] and block=3 will handle elements 266c8e769f6Sjeremylt [3*blksize : 4*blksize] 267c8e769f6Sjeremylt u: input vector 268c8e769f6Sjeremylt v: output vector 269c8e769f6Sjeremylt **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 270c8e769f6Sjeremylt 271c8e769f6Sjeremylt # libCEED call 272a8d32208Sjeremylt self._elemrestriction.apply_block(block, u, v, request=request, 273c8e769f6Sjeremylt tmode=TRANSPOSE) 274c8e769f6Sjeremylt 275c8e769f6Sjeremylt# ------------------------------------------------------------------------------ 276