1# Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors 2# All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3# 4# SPDX-License-Identifier: BSD-2-Clause 5# 6# This file is part of CEED: http://github.com/ceed 7 8from _ceed_cffi import ffi, lib 9import tempfile 10import numpy as np 11from abc import ABC 12from .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, MEM_HOST, USE_POINTER, COPY_VALUES, TRANSPOSE, NOTRANSPOSE 13from .ceed_vector import _VectorWrap 14 15# ------------------------------------------------------------------------------ 16 17 18class _ElemRestrictionBase(ABC): 19 """Ceed ElemRestriction: restriction from local vectors to elements.""" 20 21 # Destructor 22 def __del__(self): 23 # libCEED call 24 err_code = lib.CeedElemRestrictionDestroy(self._pointer) 25 self._ceed._check_error(err_code) 26 27 # Representation 28 def __repr__(self): 29 return "<CeedElemRestriction instance at " + hex(id(self)) + ">" 30 31 # String conversion for print() to stdout 32 def __str__(self): 33 """View an ElemRestriction via print().""" 34 35 # libCEED call 36 with tempfile.NamedTemporaryFile() as key_file: 37 with open(key_file.name, 'r+') as stream_file: 38 stream = ffi.cast("FILE *", stream_file) 39 40 err_code = lib.CeedElemRestrictionView(self._pointer[0], stream) 41 self._ceed._check_error(err_code) 42 43 stream_file.seek(0) 44 out_string = stream_file.read() 45 46 return out_string 47 48 # Apply CeedElemRestriction 49 def apply(self, u, v, tmode=NOTRANSPOSE, request=REQUEST_IMMEDIATE): 50 """Restrict a local vector to an element vector or apply its transpose. 51 52 Args: 53 u: input vector 54 v: output vector 55 **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE 56 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 57 58 # libCEED call 59 err_code = lib.CeedElemRestrictionApply(self._pointer[0], tmode, u._pointer[0], 60 v._pointer[0], request) 61 self._ceed._check_error(err_code) 62 63 # Transpose an ElemRestriction 64 @property 65 def T(self): 66 """Transpose an ElemRestriction.""" 67 68 return TransposeElemRestriction(self) 69 70 # Transpose an ElemRestriction 71 @property 72 def transpose(self): 73 """Transpose an ElemRestriction.""" 74 75 return TransposeElemRestriction(self) 76 77 # Create restriction vectors 78 def create_vector(self, createLvec=True, createEvec=True): 79 """Create Vectors associated with an ElemRestriction. 80 81 Args: 82 **createLvec: flag to create local vector, default True 83 **createEvec: flag to create element vector, default True 84 85 Returns: 86 [lvec, evec]: local vector and element vector, or None if flag set to false""" 87 88 # Vector pointers 89 lvecPointer = ffi.new("CeedVector *") if createLvec else ffi.NULL 90 evecPointer = ffi.new("CeedVector *") if createEvec else ffi.NULL 91 92 # libCEED call 93 err_code = lib.CeedElemRestrictionCreateVector(self._pointer[0], lvecPointer, 94 evecPointer) 95 self._ceed._check_error(err_code) 96 97 # Return vectors 98 lvec = _VectorWrap( 99 self._ceed, lvecPointer) if createLvec else None 100 evec = _VectorWrap( 101 self._ceed, evecPointer) if createEvec else None 102 103 # Return 104 return [lvec, evec] 105 106 # Get ElemRestriction multiplicity 107 def get_multiplicity(self): 108 """Get the multiplicity of nodes in an ElemRestriction. 109 110 Returns: 111 mult: local vector containing multiplicity of nodes in ElemRestriction""" 112 113 # Create mult vector 114 [mult, evec] = self.create_vector(createEvec=False) 115 mult.set_value(0) 116 117 # libCEED call 118 err_code = lib.CeedElemRestrictionGetMultiplicity( 119 self._pointer[0], mult._pointer[0]) 120 self._ceed._check_error(err_code) 121 122 # Return 123 return mult 124 125 # Get ElemRestrition Layout 126 def get_layout(self): 127 """Get the element vector layout of an ElemRestriction. 128 129 Returns: 130 layout: Vector containing layout array, stored as [nodes, components, elements]. 131 The data for node i, component j, element k in the element 132 vector is given by i*layout[0] + j*layout[1] + k*layout[2].""" 133 134 # Create output array 135 layout = np.zeros(3, dtype="int32") 136 array_pointer = ffi.cast( 137 "CeedInt (*)[3]", 138 layout.__array_interface__['data'][0]) 139 140 # libCEED call 141 err_code = lib.CeedElemRestrictionGetELayout( 142 self._pointer[0], array_pointer) 143 self._ceed._check_error(err_code) 144 145 # Return 146 return layout 147 148# ------------------------------------------------------------------------------ 149 150 151class ElemRestriction(_ElemRestrictionBase): 152 """Ceed ElemRestriction: restriction from local vectors to elements.""" 153 154 # Constructor 155 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 156 memtype=MEM_HOST, cmode=COPY_VALUES): 157 # CeedVector object 158 self._pointer = ffi.new("CeedElemRestriction *") 159 160 # Reference to Ceed 161 self._ceed = ceed 162 163 # Store array reference if needed 164 if cmode == USE_POINTER: 165 self._array_reference = offsets 166 else: 167 self._array_reference = None 168 169 # Setup the numpy array for the libCEED call 170 offsets_pointer = ffi.new("const CeedInt *") 171 offsets_pointer = ffi.cast("const CeedInt *", 172 offsets.__array_interface__['data'][0]) 173 174 # libCEED call 175 err_code = lib.CeedElemRestrictionCreate(self._ceed._pointer[0], nelem, 176 elemsize, ncomp, compstride, 177 lsize, memtype, cmode, 178 offsets_pointer, self._pointer) 179 self._ceed._check_error(err_code) 180 181# ------------------------------------------------------------------------------ 182 183 184class OrientedElemRestriction(_ElemRestrictionBase): 185 """Ceed Oriented ElemRestriction: oriented restriction from local vectors to elements.""" 186 187 # Constructor 188 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 189 orients, memtype=MEM_HOST, cmode=COPY_VALUES): 190 # CeedVector object 191 self._pointer = ffi.new("CeedElemRestriction *") 192 193 # Reference to Ceed 194 self._ceed = ceed 195 196 # Store array reference if needed 197 if cmode == USE_POINTER: 198 self._array_reference = offsets 199 self._array_reference_aux = orients 200 else: 201 self._array_reference = None 202 self._array_reference_aux = None 203 204 # Setup the numpy arrays for the libCEED call 205 offsets_pointer = ffi.new("const CeedInt *") 206 offsets_pointer = ffi.cast("const CeedInt *", 207 offsets.__array_interface__['data'][0]) 208 orients_pointer = ffi.new("const bool *") 209 orients_pointer = ffi.cast("const bool *", 210 orients.__array_interface__['data'][0]) 211 212 # libCEED call 213 err_code = lib.CeedElemRestrictionCreateOriented(self._ceed._pointer[0], nelem, 214 elemsize, ncomp, compstride, 215 lsize, memtype, cmode, 216 offsets_pointer, orients_pointer, 217 self._pointer) 218 self._ceed._check_error(err_code) 219 220# ------------------------------------------------------------------------------ 221 222 223class CurlOrientedElemRestriction(_ElemRestrictionBase): 224 """Ceed Curl Oriented ElemRestriction: curl-oriented restriction from local vectors to elements.""" 225 226 # Constructor 227 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 228 curl_orients, memtype=MEM_HOST, cmode=COPY_VALUES): 229 # CeedVector object 230 self._pointer = ffi.new("CeedElemRestriction *") 231 232 # Reference to Ceed 233 self._ceed = ceed 234 235 # Store array reference if needed 236 if cmode == USE_POINTER: 237 self._array_reference = offsets 238 self._array_reference_aux = curl_orients 239 else: 240 self._array_reference = None 241 self._array_reference_aux = None 242 243 # Setup the numpy arrays for the libCEED call 244 offsets_pointer = ffi.new("const CeedInt *") 245 offsets_pointer = ffi.cast("const CeedInt *", 246 offsets.__array_interface__['data'][0]) 247 curl_orients_pointer = ffi.new("const CeedInt8 *") 248 curl_orients_pointer = ffi.cast("const CeedInt8 *", 249 curl_orients.__array_interface__['data'][0]) 250 251 # libCEED call 252 err_code = lib.CeedElemRestrictionCreateCurlOriented(self._ceed._pointer[0], nelem, 253 elemsize, ncomp, compstride, 254 lsize, memtype, cmode, 255 offsets_pointer, 256 curl_orients_pointer, 257 self._pointer) 258 self._ceed._check_error(err_code) 259 260# ------------------------------------------------------------------------------ 261 262 263class StridedElemRestriction(_ElemRestrictionBase): 264 """Ceed Strided ElemRestriction: strided restriction from local vectors to elements.""" 265 266 # Constructor 267 def __init__(self, ceed, nelem, elemsize, ncomp, lsize, strides): 268 # CeedVector object 269 self._pointer = ffi.new("CeedElemRestriction *") 270 271 # Reference to Ceed 272 self._ceed = ceed 273 274 # Setup the numpy array for the libCEED call 275 strides_pointer = ffi.new("const CeedInt *") 276 strides_pointer = ffi.cast("const CeedInt *", 277 strides.__array_interface__['data'][0]) 278 279 # libCEED call 280 err_code = lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0], 281 nelem, elemsize, ncomp, 282 lsize, strides_pointer, 283 self._pointer) 284 self._ceed._check_error(err_code) 285 286# ------------------------------------------------------------------------------ 287 288 289class BlockedElemRestriction(_ElemRestrictionBase): 290 """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements.""" 291 292 # Constructor 293 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 294 offsets, memtype=MEM_HOST, cmode=COPY_VALUES): 295 # CeedVector object 296 self._pointer = ffi.new("CeedElemRestriction *") 297 298 # Reference to Ceed 299 self._ceed = ceed 300 301 # Setup the numpy array for the libCEED call 302 offsets_pointer = ffi.new("const CeedInt *") 303 offsets_pointer = ffi.cast("const CeedInt *", 304 offsets.__array_interface__['data'][0]) 305 306 # libCEED call 307 err_code = lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], nelem, 308 elemsize, blksize, ncomp, 309 compstride, lsize, memtype, cmode, 310 offsets_pointer, self._pointer) 311 self._ceed._check_error(err_code) 312 313 # Transpose a Blocked ElemRestriction 314 @property 315 def T(self): 316 """Transpose a Blocked ElemRestriction.""" 317 318 return TransposeBlockedElemRestriction(self) 319 320 # Transpose a Blocked ElemRestriction 321 @property 322 def transpose(self): 323 """Transpose a Blocked ElemRestriction.""" 324 325 return TransposeBlockedElemRestriction(self) 326 327 # Apply CeedElemRestriction to single block 328 def apply_block(self, block, u, v, tmode=NOTRANSPOSE, 329 request=REQUEST_IMMEDIATE): 330 """Restrict a local vector to a block of an element vector or apply its transpose. 331 332 Args: 333 block: block number to restrict to/from, i.e. block=0 will handle 334 elements [0 : blksize] and block=3 will handle elements 335 [3*blksize : 4*blksize] 336 u: input vector 337 v: output vector 338 **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE 339 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 340 341 # libCEED call 342 err_code = lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode, 343 u._pointer[0], v._pointer[0], 344 request) 345 self._ceed._check_error(err_code) 346 347# ------------------------------------------------------------------------------ 348 349 350class BlockedOrientedElemRestriction(BlockedElemRestriction): 351 """Ceed Blocked Oriented ElemRestriction: blocked oriented restriction from local vectors to elements.""" 352 353 # Constructor 354 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 355 offsets, orients, memtype=MEM_HOST, cmode=COPY_VALUES): 356 # CeedVector object 357 self._pointer = ffi.new("CeedElemRestriction *") 358 359 # Reference to Ceed 360 self._ceed = ceed 361 362 # Setup the numpy array for the libCEED call 363 offsets_pointer = ffi.new("const CeedInt *") 364 offsets_pointer = ffi.cast("const CeedInt *", 365 offsets.__array_interface__['data'][0]) 366 orients_pointer = ffi.new("const bool *") 367 orients_pointer = ffi.cast("const bool *", 368 orients.__array_interface__['data'][0]) 369 370 # libCEED call 371 err_code = lib.CeedElemRestrictionCreateBlockedOriented(self._ceed._pointer[0], nelem, 372 elemsize, blksize, ncomp, 373 compstride, lsize, memtype, cmode, 374 offsets_pointer, orients_pointer, 375 self._pointer) 376 self._ceed._check_error(err_code) 377 378# ------------------------------------------------------------------------------ 379 380 381class BlockedCurlOrientedElemRestriction(BlockedElemRestriction): 382 """Ceed Blocked Curl Oriented ElemRestriction: blocked curl-oriented restriction from local vectors to elements.""" 383 384 # Constructor 385 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 386 offsets, curl_orients, memtype=MEM_HOST, cmode=COPY_VALUES): 387 # CeedVector object 388 self._pointer = ffi.new("CeedElemRestriction *") 389 390 # Reference to Ceed 391 self._ceed = ceed 392 393 # Setup the numpy array for the libCEED call 394 offsets_pointer = ffi.new("const CeedInt *") 395 offsets_pointer = ffi.cast("const CeedInt *", 396 offsets.__array_interface__['data'][0]) 397 curl_orients_pointer = ffi.new("const CeedInt8 *") 398 curl_orients_pointer = ffi.cast("const CeedInt8 *", 399 curl_orients.__array_interface__['data'][0]) 400 401 # libCEED call 402 err_code = lib.CeedElemRestrictionCreateBlockedCurlOriented(self._ceed._pointer[0], nelem, 403 elemsize, blksize, ncomp, 404 compstride, lsize, memtype, cmode, 405 offsets_pointer, curl_orients_pointer, 406 self._pointer) 407 self._ceed._check_error(err_code) 408 409# ------------------------------------------------------------------------------ 410 411 412class BlockedStridedElemRestriction(BlockedElemRestriction): 413 """Ceed Blocked Strided ElemRestriction: blocked strided restriction from local vectors to elements.""" 414 415 # Constructor 416 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, lsize, strides): 417 # CeedVector object 418 self._pointer = ffi.new("CeedElemRestriction *") 419 420 # Reference to Ceed 421 self._ceed = ceed 422 423 # Setup the numpy array for the libCEED call 424 strides_pointer = ffi.new("const CeedInt *") 425 strides_pointer = ffi.cast("const CeedInt *", 426 strides.__array_interface__['data'][0]) 427 428 # libCEED call 429 err_code = lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem, 430 elemsize, blksize, ncomp, 431 lsize, strides_pointer, 432 self._pointer) 433 self._ceed._check_error(err_code) 434 435# ------------------------------------------------------------------------------ 436 437 438class TransposeElemRestriction(): 439 """Ceed ElemRestriction: transpose restriction from elements to local vectors.""" 440 441 # Attributes 442 _elemrestriction = None 443 444 # Constructor 445 def __init__(self, elemrestriction): 446 447 # Reference elemrestriction 448 self._elemrestriction = elemrestriction 449 450 # Representation 451 def __repr__(self): 452 return "<Transpose CeedElemRestriction instance at " + \ 453 hex(id(self)) + ">" 454 455 # Apply Transpose CeedElemRestriction 456 457 def apply(self, u, v, request=REQUEST_IMMEDIATE): 458 """Restrict an element vector to a local vector. 459 460 Args: 461 u: input vector 462 v: output vector 463 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 464 465 # libCEED call 466 self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE) 467 468# ------------------------------------------------------------------------------ 469 470 471class TransposeBlockedElemRestriction(TransposeElemRestriction): 472 """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements 473 to local vectors.""" 474 475 # Apply Transpose CeedElemRestriction 476 def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE): 477 """Restrict a block of an element vector to a local vector. 478 479 Args: 480 block: block number to restrict to/from, i.e. block=0 will handle 481 elements [0 : blksize] and block=3 will handle elements 482 [3*blksize : 4*blksize] 483 u: input vector 484 v: output vector 485 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 486 487 # libCEED call 488 self._elemrestriction.apply_block(block, u, v, request=request, 489 tmode=TRANSPOSE) 490 491# ------------------------------------------------------------------------------ 492