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 L-vector Layout 126 def get_l_layout(self): 127 """Get the local 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 layout_pointer = ffi.cast("const CeedInt *", 137 layout.__array_interface__['data'][0]) 138 139 # libCEED call 140 err_code = lib.CeedElemRestrictionGetLLayout( 141 self._pointer[0], layout_pointer) 142 self._ceed._check_error(err_code) 143 144 # Return 145 return layout 146 147 # Get ElemRestrition E-vector Layout 148 def get_e_layout(self): 149 """Get the element vector layout of an ElemRestriction. 150 151 Returns: 152 layout: Vector containing layout array, stored as [nodes, components, elements]. 153 The data for node i, component j, element k in the element 154 vector is given by i*layout[0] + j*layout[1] + k*layout[2].""" 155 156 # Create output array 157 layout = np.zeros(3, dtype="int32") 158 layout_pointer = ffi.cast("const CeedInt *", 159 layout.__array_interface__['data'][0]) 160 161 # libCEED call 162 err_code = lib.CeedElemRestrictionGetELayout( 163 self._pointer[0], layout_pointer) 164 self._ceed._check_error(err_code) 165 166 # Return 167 return layout 168 169# ------------------------------------------------------------------------------ 170 171 172class ElemRestriction(_ElemRestrictionBase): 173 """Ceed ElemRestriction: restriction from local vectors to elements.""" 174 175 # Constructor 176 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 177 memtype=MEM_HOST, cmode=COPY_VALUES): 178 # CeedVector object 179 self._pointer = ffi.new("CeedElemRestriction *") 180 181 # Reference to Ceed 182 self._ceed = ceed 183 184 # Store array reference if needed 185 if cmode == USE_POINTER: 186 self._array_reference = offsets 187 else: 188 self._array_reference = None 189 190 # Setup the numpy array for the libCEED call 191 offsets_pointer = ffi.new("const CeedInt *") 192 offsets_pointer = ffi.cast("const CeedInt *", 193 offsets.__array_interface__['data'][0]) 194 195 # libCEED call 196 err_code = lib.CeedElemRestrictionCreate(self._ceed._pointer[0], nelem, 197 elemsize, ncomp, compstride, 198 lsize, memtype, cmode, 199 offsets_pointer, self._pointer) 200 self._ceed._check_error(err_code) 201 202# ------------------------------------------------------------------------------ 203 204 205class OrientedElemRestriction(_ElemRestrictionBase): 206 """Ceed Oriented ElemRestriction: oriented restriction from local vectors to elements.""" 207 208 # Constructor 209 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 210 orients, memtype=MEM_HOST, cmode=COPY_VALUES): 211 # CeedVector object 212 self._pointer = ffi.new("CeedElemRestriction *") 213 214 # Reference to Ceed 215 self._ceed = ceed 216 217 # Store array reference if needed 218 if cmode == USE_POINTER: 219 self._array_reference = offsets 220 self._array_reference_aux = orients 221 else: 222 self._array_reference = None 223 self._array_reference_aux = None 224 225 # Setup the numpy arrays for the libCEED call 226 offsets_pointer = ffi.new("const CeedInt *") 227 offsets_pointer = ffi.cast("const CeedInt *", 228 offsets.__array_interface__['data'][0]) 229 orients_pointer = ffi.new("const bool *") 230 orients_pointer = ffi.cast("const bool *", 231 orients.__array_interface__['data'][0]) 232 233 # libCEED call 234 err_code = lib.CeedElemRestrictionCreateOriented(self._ceed._pointer[0], nelem, 235 elemsize, ncomp, compstride, 236 lsize, memtype, cmode, 237 offsets_pointer, orients_pointer, 238 self._pointer) 239 self._ceed._check_error(err_code) 240 241# ------------------------------------------------------------------------------ 242 243 244class CurlOrientedElemRestriction(_ElemRestrictionBase): 245 """Ceed Curl Oriented ElemRestriction: curl-oriented restriction from local vectors to elements.""" 246 247 # Constructor 248 def __init__(self, ceed, nelem, elemsize, ncomp, compstride, lsize, offsets, 249 curl_orients, memtype=MEM_HOST, cmode=COPY_VALUES): 250 # CeedVector object 251 self._pointer = ffi.new("CeedElemRestriction *") 252 253 # Reference to Ceed 254 self._ceed = ceed 255 256 # Store array reference if needed 257 if cmode == USE_POINTER: 258 self._array_reference = offsets 259 self._array_reference_aux = curl_orients 260 else: 261 self._array_reference = None 262 self._array_reference_aux = None 263 264 # Setup the numpy arrays for the libCEED call 265 offsets_pointer = ffi.new("const CeedInt *") 266 offsets_pointer = ffi.cast("const CeedInt *", 267 offsets.__array_interface__['data'][0]) 268 curl_orients_pointer = ffi.new("const CeedInt8 *") 269 curl_orients_pointer = ffi.cast("const CeedInt8 *", 270 curl_orients.__array_interface__['data'][0]) 271 272 # libCEED call 273 err_code = lib.CeedElemRestrictionCreateCurlOriented(self._ceed._pointer[0], nelem, 274 elemsize, ncomp, compstride, 275 lsize, memtype, cmode, 276 offsets_pointer, 277 curl_orients_pointer, 278 self._pointer) 279 self._ceed._check_error(err_code) 280 281# ------------------------------------------------------------------------------ 282 283 284class StridedElemRestriction(_ElemRestrictionBase): 285 """Ceed Strided ElemRestriction: strided restriction from local vectors to elements.""" 286 287 # Constructor 288 def __init__(self, ceed, nelem, elemsize, ncomp, lsize, strides): 289 # CeedVector object 290 self._pointer = ffi.new("CeedElemRestriction *") 291 292 # Reference to Ceed 293 self._ceed = ceed 294 295 # Setup the numpy array for the libCEED call 296 strides_pointer = ffi.new("const CeedInt *") 297 strides_pointer = ffi.cast("const CeedInt *", 298 strides.__array_interface__['data'][0]) 299 300 # libCEED call 301 err_code = lib.CeedElemRestrictionCreateStrided(self._ceed._pointer[0], 302 nelem, elemsize, ncomp, 303 lsize, strides_pointer, 304 self._pointer) 305 self._ceed._check_error(err_code) 306 307# ------------------------------------------------------------------------------ 308 309 310class BlockedElemRestriction(_ElemRestrictionBase): 311 """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements.""" 312 313 # Constructor 314 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 315 offsets, memtype=MEM_HOST, cmode=COPY_VALUES): 316 # CeedVector object 317 self._pointer = ffi.new("CeedElemRestriction *") 318 319 # Reference to Ceed 320 self._ceed = ceed 321 322 # Setup the numpy array for the libCEED call 323 offsets_pointer = ffi.new("const CeedInt *") 324 offsets_pointer = ffi.cast("const CeedInt *", 325 offsets.__array_interface__['data'][0]) 326 327 # libCEED call 328 err_code = lib.CeedElemRestrictionCreateBlocked(self._ceed._pointer[0], nelem, 329 elemsize, blksize, ncomp, 330 compstride, lsize, memtype, cmode, 331 offsets_pointer, self._pointer) 332 self._ceed._check_error(err_code) 333 334 # Transpose a Blocked ElemRestriction 335 @property 336 def T(self): 337 """Transpose a Blocked ElemRestriction.""" 338 339 return TransposeBlockedElemRestriction(self) 340 341 # Transpose a Blocked ElemRestriction 342 @property 343 def transpose(self): 344 """Transpose a Blocked ElemRestriction.""" 345 346 return TransposeBlockedElemRestriction(self) 347 348 # Apply CeedElemRestriction to single block 349 def apply_block(self, block, u, v, tmode=NOTRANSPOSE, 350 request=REQUEST_IMMEDIATE): 351 """Restrict a local vector to a block of an element vector or apply its transpose. 352 353 Args: 354 block: block number to restrict to/from, i.e. block=0 will handle 355 elements [0 : blksize] and block=3 will handle elements 356 [3*blksize : 4*blksize] 357 u: input vector 358 v: output vector 359 **tmode: apply restriction or transpose, default CEED_NOTRANSPOSE 360 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 361 362 # libCEED call 363 err_code = lib.CeedElemRestrictionApplyBlock(self._pointer[0], block, tmode, 364 u._pointer[0], v._pointer[0], 365 request) 366 self._ceed._check_error(err_code) 367 368# ------------------------------------------------------------------------------ 369 370 371class BlockedOrientedElemRestriction(BlockedElemRestriction): 372 """Ceed Blocked Oriented ElemRestriction: blocked oriented restriction from local vectors to elements.""" 373 374 # Constructor 375 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 376 offsets, orients, memtype=MEM_HOST, cmode=COPY_VALUES): 377 # CeedVector object 378 self._pointer = ffi.new("CeedElemRestriction *") 379 380 # Reference to Ceed 381 self._ceed = ceed 382 383 # Setup the numpy array for the libCEED call 384 offsets_pointer = ffi.new("const CeedInt *") 385 offsets_pointer = ffi.cast("const CeedInt *", 386 offsets.__array_interface__['data'][0]) 387 orients_pointer = ffi.new("const bool *") 388 orients_pointer = ffi.cast("const bool *", 389 orients.__array_interface__['data'][0]) 390 391 # libCEED call 392 err_code = lib.CeedElemRestrictionCreateBlockedOriented(self._ceed._pointer[0], nelem, 393 elemsize, blksize, ncomp, 394 compstride, lsize, memtype, cmode, 395 offsets_pointer, orients_pointer, 396 self._pointer) 397 self._ceed._check_error(err_code) 398 399# ------------------------------------------------------------------------------ 400 401 402class BlockedCurlOrientedElemRestriction(BlockedElemRestriction): 403 """Ceed Blocked Curl Oriented ElemRestriction: blocked curl-oriented restriction from local vectors to elements.""" 404 405 # Constructor 406 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, compstride, lsize, 407 offsets, curl_orients, memtype=MEM_HOST, cmode=COPY_VALUES): 408 # CeedVector object 409 self._pointer = ffi.new("CeedElemRestriction *") 410 411 # Reference to Ceed 412 self._ceed = ceed 413 414 # Setup the numpy array for the libCEED call 415 offsets_pointer = ffi.new("const CeedInt *") 416 offsets_pointer = ffi.cast("const CeedInt *", 417 offsets.__array_interface__['data'][0]) 418 curl_orients_pointer = ffi.new("const CeedInt8 *") 419 curl_orients_pointer = ffi.cast("const CeedInt8 *", 420 curl_orients.__array_interface__['data'][0]) 421 422 # libCEED call 423 err_code = lib.CeedElemRestrictionCreateBlockedCurlOriented(self._ceed._pointer[0], nelem, 424 elemsize, blksize, ncomp, 425 compstride, lsize, memtype, cmode, 426 offsets_pointer, curl_orients_pointer, 427 self._pointer) 428 self._ceed._check_error(err_code) 429 430# ------------------------------------------------------------------------------ 431 432 433class BlockedStridedElemRestriction(BlockedElemRestriction): 434 """Ceed Blocked Strided ElemRestriction: blocked strided restriction from local vectors to elements.""" 435 436 # Constructor 437 def __init__(self, ceed, nelem, elemsize, blksize, ncomp, lsize, strides): 438 # CeedVector object 439 self._pointer = ffi.new("CeedElemRestriction *") 440 441 # Reference to Ceed 442 self._ceed = ceed 443 444 # Setup the numpy array for the libCEED call 445 strides_pointer = ffi.new("const CeedInt *") 446 strides_pointer = ffi.cast("const CeedInt *", 447 strides.__array_interface__['data'][0]) 448 449 # libCEED call 450 err_code = lib.CeedElemRestrictionCreateBlockedStrided(self._ceed._pointer[0], nelem, 451 elemsize, blksize, ncomp, 452 lsize, strides_pointer, 453 self._pointer) 454 self._ceed._check_error(err_code) 455 456# ------------------------------------------------------------------------------ 457 458 459class TransposeElemRestriction(): 460 """Ceed ElemRestriction: transpose restriction from elements to local vectors.""" 461 462 # Attributes 463 _elemrestriction = None 464 465 # Constructor 466 def __init__(self, elemrestriction): 467 468 # Reference elemrestriction 469 self._elemrestriction = elemrestriction 470 471 # Representation 472 def __repr__(self): 473 return "<Transpose CeedElemRestriction instance at " + \ 474 hex(id(self)) + ">" 475 476 # Apply Transpose CeedElemRestriction 477 478 def apply(self, u, v, request=REQUEST_IMMEDIATE): 479 """Restrict an element vector to a local vector. 480 481 Args: 482 u: input vector 483 v: output vector 484 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 485 486 # libCEED call 487 self._elemrestriction.apply(u, v, request=request, tmode=TRANSPOSE) 488 489# ------------------------------------------------------------------------------ 490 491 492class TransposeBlockedElemRestriction(TransposeElemRestriction): 493 """Transpose Ceed Blocked ElemRestriction: blocked transpose restriction from elements 494 to local vectors.""" 495 496 # Apply Transpose CeedElemRestriction 497 def apply_block(self, block, u, v, request=REQUEST_IMMEDIATE): 498 """Restrict a block of an element vector to a local vector. 499 500 Args: 501 block: block number to restrict to/from, i.e. block=0 will handle 502 elements [0 : blksize] and block=3 will handle elements 503 [3*blksize : 4*blksize] 504 u: input vector 505 v: output vector 506 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 507 508 # libCEED call 509 self._elemrestriction.apply_block(block, u, v, request=request, 510 tmode=TRANSPOSE) 511 512# ------------------------------------------------------------------------------ 513