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 sys 10import os 11import io 12import numpy as np 13import tempfile 14from abc import ABC 15from .ceed_vector import Vector 16from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1, BasisHdiv, BasisHcurl 17from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction 18from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction 19from .ceed_qfunctioncontext import QFunctionContext 20from .ceed_operator import Operator, CompositeOperator 21from .ceed_constants import * 22 23# ------------------------------------------------------------------------------ 24 25 26class Ceed(): 27 """Ceed: core components.""" 28 29 # Constructor 30 def __init__(self, resource="/cpu/self", on_error="store"): 31 # libCEED object 32 self._pointer = ffi.new("Ceed *") 33 34 # libCEED call 35 resourceAscii = ffi.new("char[]", resource.encode("ascii")) 36 os.environ["CEED_ERROR_HANDLER"] = "return" 37 err_code = lib.CeedInit(resourceAscii, self._pointer) 38 if err_code: 39 raise Exception("Error initializing backend resource: " + resource) 40 error_handlers = dict( 41 store="CeedErrorStore", 42 abort="CeedErrorAbort", 43 ) 44 lib.CeedSetErrorHandler( 45 self._pointer[0], ffi.addressof( 46 lib, error_handlers[on_error])) 47 48 # Representation 49 def __repr__(self): 50 return "<Ceed instance at " + hex(id(self)) + ">" 51 52 # String conversion for print() to stdout 53 def __str__(self): 54 """View a Ceed via print().""" 55 56 # libCEED call 57 with tempfile.NamedTemporaryFile() as key_file: 58 with open(key_file.name, 'r+') as stream_file: 59 stream = ffi.cast("FILE *", stream_file) 60 61 err_code = lib.CeedView(self._pointer[0], stream) 62 self._check_error(err_code) 63 64 stream_file.seek(0) 65 out_string = stream_file.read() 66 67 return out_string 68 69 # Error handler 70 def _check_error(self, err_code): 71 """Check return code and retrieve error message for non-zero code""" 72 if (err_code != lib.CEED_ERROR_SUCCESS): 73 message = ffi.new("char **") 74 lib.CeedGetErrorMessage(self._pointer[0], message) 75 raise Exception(ffi.string(message[0]).decode("UTF-8")) 76 77 # Get Resource 78 def get_resource(self): 79 """Get the full resource name for a Ceed context. 80 81 Returns: 82 resource: resource name""" 83 84 # libCEED call 85 resource = ffi.new("char **") 86 err_code = lib.CeedGetResource(self._pointer[0], resource) 87 self._check_error(err_code) 88 89 return ffi.string(resource[0]).decode("UTF-8") 90 91 # Preferred MemType 92 def get_preferred_memtype(self): 93 """Return Ceed preferred memory type. 94 95 Returns: 96 memtype: Ceed preferred memory type""" 97 98 # libCEED call 99 memtype = ffi.new("CeedMemType *", MEM_HOST) 100 err_code = lib.CeedGetPreferredMemType(self._pointer[0], memtype) 101 self._check_error(err_code) 102 103 return memtype[0] 104 105 # Convenience function to get CeedScalar type 106 def scalar_type(self): 107 """Return dtype string for CeedScalar. 108 109 Returns: 110 np_dtype: String naming numpy data type corresponding to CeedScalar""" 111 112 return scalar_types[lib.CEED_SCALAR_TYPE] 113 114 # --- Basis utility functions --- 115 116 # Gauss quadrature 117 def gauss_quadrature(self, q): 118 """Construct a Gauss-Legendre quadrature. 119 120 Args: 121 Q: number of quadrature points (integrates polynomials of 122 degree 2*Q-1 exactly) 123 124 Returns: 125 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 126 and array of length Q to hold the weights""" 127 128 # Setup arguments 129 qref1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 130 qweight1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 131 132 qref1d_pointer = ffi.new("CeedScalar *") 133 qref1d_pointer = ffi.cast( 134 "CeedScalar *", 135 qref1d.__array_interface__['data'][0]) 136 137 qweight1d_pointer = ffi.new("CeedScalar *") 138 qweight1d_pointer = ffi.cast( 139 "CeedScalar *", 140 qweight1d.__array_interface__['data'][0]) 141 142 # libCEED call 143 err_code = lib.CeedGaussQuadrature(q, qref1d_pointer, qweight1d_pointer) 144 self._check_error(err_code) 145 146 return qref1d, qweight1d 147 148 # Lobatto quadrature 149 def lobatto_quadrature(self, q): 150 """Construct a Gauss-Legendre-Lobatto quadrature. 151 152 Args: 153 q: number of quadrature points (integrates polynomials of 154 degree 2*Q-3 exactly) 155 156 Returns: 157 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 158 and array of length Q to hold the weights""" 159 160 # Setup arguments 161 qref1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 162 qref1d_pointer = ffi.new("CeedScalar *") 163 qref1d_pointer = ffi.cast( 164 "CeedScalar *", 165 qref1d.__array_interface__['data'][0]) 166 167 qweight1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 168 qweight1d_pointer = ffi.new("CeedScalar *") 169 qweight1d_pointer = ffi.cast( 170 "CeedScalar *", 171 qweight1d.__array_interface__['data'][0]) 172 173 # libCEED call 174 err_code = lib.CeedLobattoQuadrature( 175 q, qref1d_pointer, qweight1d_pointer) 176 self._check_error(err_code) 177 178 return qref1d, qweight1d 179 180 # --- libCEED objects --- 181 182 # CeedVector 183 def Vector(self, size): 184 """Ceed Vector: storing and manipulating vectors. 185 186 Args: 187 size: length of vector 188 189 Returns: 190 vector: Ceed Vector""" 191 192 return Vector(self, size) 193 194 # CeedElemRestriction 195 def ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, offsets, 196 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES): 197 """Ceed ElemRestriction: restriction from local vectors to elements. 198 199 Args: 200 nelem: number of elements described by the restriction 201 elemsize: size (number of nodes) per element 202 ncomp: number of field components per interpolation node 203 (1 for scalar fields) 204 compstride: Stride between components for the same L-vector "node". 205 Data for node i, component k can be found in the 206 L-vector at index [offsets[i] + k*compstride]. 207 lsize: The size of the L-vector. This vector may be larger than 208 the elements and fields given by this restriction. 209 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 210 ordered list of the offsets (into the input Ceed Vector) 211 for the unknowns corresponding to element i, where 212 0 <= i < nelem. All offsets must be in the range 213 [0, lsize - 1]. 214 **memtype: memory type of the offsets array, default CEED_MEM_HOST 215 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 216 217 Returns: 218 elemrestriction: Ceed ElemRestiction""" 219 220 return ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, 221 offsets, memtype=memtype, cmode=cmode) 222 223 def StridedElemRestriction(self, nelem, elemsize, ncomp, lsize, strides): 224 """Ceed Identity ElemRestriction: strided restriction from local vectors 225 to elements. 226 227 Args: 228 nelem: number of elements described by the restriction 229 elemsize: size (number of nodes) per element 230 ncomp: number of field components per interpolation node 231 (1 for scalar fields) 232 lsize: The size of the L-vector. This vector may be larger than 233 the elements and fields given by this restriction. 234 *strides: Array for strides between [nodes, components, elements]. 235 The data for node i, component j, element k in the 236 L-vector is given by 237 i*strides[0] + j*strides[1] + k*strides[2] 238 239 Returns: 240 elemrestriction: Ceed Strided ElemRestiction""" 241 242 return StridedElemRestriction( 243 self, nelem, elemsize, ncomp, lsize, strides) 244 245 def BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, compstride, 246 lsize, offsets, memtype=lib.CEED_MEM_HOST, 247 cmode=lib.CEED_COPY_VALUES): 248 """Ceed Blocked ElemRestriction: blocked restriction from local vectors 249 to elements. 250 251 Args: 252 nelem: number of elements described by the restriction 253 elemsize: size (number of nodes) per element 254 blksize: number of elements in a block 255 ncomp: number of field components per interpolation node 256 (1 for scalar fields) 257 lsize: The size of the L-vector. This vector may be larger than 258 the elements and fields given by this restriction. 259 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 260 ordered list of the offsets (into the input Ceed Vector) 261 for the unknowns corresponding to element i, where 262 0 <= i < nelem. All offsets must be in the range 263 [0, lsize - 1]. The backend will permute and pad this 264 array to the desired ordering for the blocksize, which is 265 typically given by the backend. The default reordering is 266 to interlace elements. 267 **memtype: memory type of the offsets array, default CEED_MEM_HOST 268 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 269 270 Returns: 271 elemrestriction: Ceed Blocked ElemRestiction""" 272 273 return BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, 274 compstride, lsize, offsets, 275 memtype=memtype, cmode=cmode) 276 277 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, ncomp, 278 lsize, strides): 279 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction 280 from local vectors to elements. 281 282 Args: 283 nelem: number of elements described in the restriction 284 elemsize: size (number of nodes) per element 285 blksize: number of elements in a block 286 ncomp: number of field components per interpolation node 287 (1 for scalar fields) 288 lsize: The size of the L-vector. This vector may be larger than 289 the elements and fields given by this restriction. 290 *strides: Array for strides between [nodes, components, elements]. 291 The data for node i, component j, element k in the 292 L-vector is given by 293 i*strides[0] + j*strides[1] + k*strides[2] 294 295 Returns: 296 elemrestriction: Ceed Strided ElemRestiction""" 297 298 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 299 ncomp, lsize, strides) 300 301 # CeedBasis 302 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 303 qref1d, qweight1d): 304 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 305 H^1 discretizations. 306 307 Args: 308 dim: topological dimension 309 ncomp: number of field components (1 for scalar fields) 310 P1d: number of nodes in one dimension 311 Q1d: number of quadrature points in one dimension 312 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix 313 expressing the values of nodal basis functions at 314 quadrature points 315 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix 316 expressing the derivatives of nodal basis functions at 317 quadrature points 318 *qref1d: Numpy array of length Q1d holding the locations of 319 quadrature points on the 1D reference element [-1, 1] 320 *qweight1d: Numpy array of length Q1d holding the quadrature 321 weights on the reference element 322 323 Returns: 324 basis: Ceed Basis""" 325 326 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 327 qref1d, qweight1d) 328 329 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 330 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange 331 basis objects for H^1 discretizations. 332 333 Args: 334 dim: topological dimension 335 ncomp: number of field components (1 for scalar fields) 336 P: number of Gauss-Lobatto nodes in one dimension. The 337 polynomial degree of the resulting Q_k element is k=P-1. 338 Q: number of quadrature points in one dimension 339 qmode: distribution of the Q quadrature points (affects order of 340 accuracy for the quadrature) 341 342 Returns: 343 basis: Ceed Basis""" 344 345 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 346 347 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 348 """Ceed H1 Basis: finite element non tensor-product basis for H^1 349 discretizations. 350 351 Args: 352 topo: topology of the element, e.g. hypercube, simplex, etc 353 ncomp: number of field components (1 for scalar fields) 354 nnodes: total number of nodes 355 nqpts: total number of quadrature points 356 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix 357 expressing the values of nodal basis functions at 358 quadrature points 359 *grad: Numpy array holding the row-major (dim * nqpts * nnodes) 360 matrix expressing the derivatives of nodal basis functions 361 at quadrature points 362 *qref: Numpy array of length (nqpts * dim) holding the locations of 363 quadrature points on the reference element [-1, 1] 364 *qweight: Numpy array of length nnodes holding the quadrature 365 weights on the reference element 366 367 Returns: 368 basis: Ceed Basis""" 369 370 return BasisH1(self, topo, ncomp, nnodes, nqpts, 371 interp, grad, qref, qweight) 372 373 def BasisHdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight): 374 """Ceed Hdiv Basis: finite element non tensor-product basis for H(div) 375 discretizations. 376 377 Args: 378 topo: topology of the element, e.g. hypercube, simplex, etc 379 ncomp: number of field components (1 for scalar fields) 380 nnodes: total number of nodes 381 nqpts: total number of quadrature points 382 *interp: Numpy array holding the row-major (dim * nqpts * nnodes) 383 matrix expressing the values of basis functions at 384 quadrature points 385 *div: Numpy array holding the row-major (nqpts * nnodes) matrix 386 expressing the divergence of basis functions at 387 quadrature points 388 *qref: Numpy array of length (nqpts * dim) holding the locations of 389 quadrature points on the reference element [-1, 1] 390 *qweight: Numpy array of length nnodes holding the quadrature 391 weights on the reference element 392 393 Returns: 394 basis: Ceed Basis""" 395 396 return BasisHdiv(self, topo, ncomp, nnodes, nqpts, 397 interp, div, qref, qweight) 398 399 def BasisHcurl(self, topo, ncomp, nnodes, nqpts, 400 interp, curl, qref, qweight): 401 """Ceed Hcurl Basis: finite element non tensor-product basis for H(curl) 402 discretizations. 403 404 Args: 405 topo: topology of the element, e.g. hypercube, simplex, etc 406 ncomp: number of field components (1 for scalar fields) 407 nnodes: total number of nodes 408 nqpts: total number of quadrature points 409 *interp: Numpy array holding the row-major (dim * nqpts * nnodes) 410 matrix expressing the values of basis functions at 411 quadrature points 412 *curl: Numpy array holding the row-major (curlcomp * nqpts * nnodes), 413 curlcomp = 1 if dim < 3 else dim, matrix expressing the curl 414 of basis functions at quadrature points 415 *qref: Numpy array of length (nqpts * dim) holding the locations of 416 quadrature points on the reference element [-1, 1] 417 *qweight: Numpy array of length nnodes holding the quadrature 418 weights on the reference element 419 420 Returns: 421 basis: Ceed Basis""" 422 423 return BasisHcurl(self, topo, ncomp, nnodes, nqpts, 424 interp, curl, qref, qweight) 425 426 # CeedQFunction 427 def QFunction(self, vlength, f, source): 428 """Ceed QFunction: point-wise operation at quadrature points for 429 evaluating volumetric terms. 430 431 Args: 432 vlength: vector length. Caller must ensure that number of quadrature 433 points is a multiple of vlength 434 f: ctypes function pointer to evaluate action at quadrature points 435 source: absolute path to source of QFunction, 436 "\\abs_path\\file.h:function_name 437 438 Returns: 439 qfunction: Ceed QFunction""" 440 441 return QFunction(self, vlength, f, source) 442 443 def QFunctionByName(self, name): 444 """Ceed QFunction By Name: point-wise operation at quadrature points 445 from a given gallery, for evaluating volumetric terms. 446 447 Args: 448 name: name of QFunction to use from gallery 449 450 Returns: 451 qfunction: Ceed QFunction By Name""" 452 453 return QFunctionByName(self, name) 454 455 def IdentityQFunction(self, size, inmode, outmode): 456 """Ceed Idenity QFunction: identity qfunction operation. 457 458 Args: 459 size: size of the qfunction fields 460 **inmode: CeedEvalMode for input to Ceed QFunction 461 **outmode: CeedEvalMode for output to Ceed QFunction 462 463 Returns: 464 qfunction: Ceed Identity QFunction""" 465 466 return IdentityQFunction(self, size, inmode, outmode) 467 468 def QFunctionContext(self): 469 """Ceed QFunction Context: stores Ceed QFunction user context data. 470 471 Returns: 472 userContext: Ceed QFunction Context""" 473 474 return QFunctionContext(self) 475 476 # CeedOperator 477 def Operator(self, qf, dqf=None, qdfT=None): 478 """Ceed Operator: composed FE-type operations on vectors. 479 480 Args: 481 qf: QFunction defining the action of the operator at quadrature 482 points 483 **dqf: QFunction defining the action of the Jacobian, default None 484 **dqfT: QFunction defining the action of the transpose of the 485 Jacobian, default None 486 487 Returns: 488 operator: Ceed Operator""" 489 490 return Operator(self, qf, dqf, qdfT) 491 492 def CompositeOperator(self): 493 """Composite Ceed Operator: composition of multiple CeedOperators. 494 495 Returns: 496 operator: Ceed Composite Operator""" 497 498 return CompositeOperator(self) 499 500 # Destructor 501 def __del__(self): 502 # libCEED call 503 lib.CeedDestroy(self._pointer) 504 505# ------------------------------------------------------------------------------ 506