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 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 (nqpts * dim * 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 # CeedQFunction 374 def QFunction(self, vlength, f, source): 375 """Ceed QFunction: point-wise operation at quadrature points for 376 evaluating volumetric terms. 377 378 Args: 379 vlength: vector length. Caller must ensure that number of quadrature 380 points is a multiple of vlength 381 f: ctypes function pointer to evaluate action at quadrature points 382 source: absolute path to source of QFunction, 383 "\\abs_path\\file.h:function_name 384 385 Returns: 386 qfunction: Ceed QFunction""" 387 388 return QFunction(self, vlength, f, source) 389 390 def QFunctionByName(self, name): 391 """Ceed QFunction By Name: point-wise operation at quadrature points 392 from a given gallery, for evaluating volumetric terms. 393 394 Args: 395 name: name of QFunction to use from gallery 396 397 Returns: 398 qfunction: Ceed QFunction By Name""" 399 400 return QFunctionByName(self, name) 401 402 def IdentityQFunction(self, size, inmode, outmode): 403 """Ceed Idenity QFunction: identity qfunction operation. 404 405 Args: 406 size: size of the qfunction fields 407 **inmode: CeedEvalMode for input to Ceed QFunction 408 **outmode: CeedEvalMode for output to Ceed QFunction 409 410 Returns: 411 qfunction: Ceed Identity QFunction""" 412 413 return IdentityQFunction(self, size, inmode, outmode) 414 415 def QFunctionContext(self): 416 """Ceed QFunction Context: stores Ceed QFunction user context data. 417 418 Returns: 419 userContext: Ceed QFunction Context""" 420 421 return QFunctionContext(self) 422 423 # CeedOperator 424 def Operator(self, qf, dqf=None, qdfT=None): 425 """Ceed Operator: composed FE-type operations on vectors. 426 427 Args: 428 qf: QFunction defining the action of the operator at quadrature 429 points 430 **dqf: QFunction defining the action of the Jacobian, default None 431 **dqfT: QFunction defining the action of the transpose of the 432 Jacobian, default None 433 434 Returns: 435 operator: Ceed Operator""" 436 437 return Operator(self, qf, dqf, qdfT) 438 439 def CompositeOperator(self): 440 """Composite Ceed Operator: composition of multiple CeedOperators. 441 442 Returns: 443 operator: Ceed Composite Operator""" 444 445 return CompositeOperator(self) 446 447 # Destructor 448 def __del__(self): 449 # libCEED call 450 lib.CeedDestroy(self._pointer) 451 452# ------------------------------------------------------------------------------ 453