1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17from _ceed_cffi import ffi, lib 18import sys 19import io 20import tempfile 21from abc import ABC 22from .ceed_vector import Vector 23from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1 24from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction 25from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction 26from .ceed_operator import Operator, CompositeOperator 27from .ceed_constants import * 28 29# ------------------------------------------------------------------------------ 30 31 32class Ceed(): 33 """Ceed: core components.""" 34 # Attributes 35 _pointer = ffi.NULL 36 37 # Constructor 38 def __init__(self, resource="/cpu/self"): 39 # libCEED object 40 self._pointer = ffi.new("Ceed *") 41 42 # libCEED call 43 resourceAscii = ffi.new("char[]", resource.encode("ascii")) 44 lib.CeedInit(resourceAscii, self._pointer) 45 46 # Representation 47 def __repr__(self): 48 return "<Ceed instance at " + hex(id(self)) + ">" 49 50 # String conversion for print() to stdout 51 def __str__(self): 52 """View a Ceed via print().""" 53 54 # libCEED call 55 with tempfile.NamedTemporaryFile() as key_file: 56 with open(key_file.name, 'r+') as stream_file: 57 stream = ffi.cast("FILE *", stream_file) 58 59 lib.CeedView(self._pointer[0], stream) 60 61 stream_file.seek(0) 62 out_string = stream_file.read() 63 64 return out_string 65 66 # Get Resource 67 def get_resource(self): 68 """Get the full resource name for a Ceed context. 69 70 Returns: 71 resource: resource name""" 72 73 # libCEED call 74 resource = ffi.new("char **") 75 lib.CeedGetResource(self._pointer[0], resource) 76 77 return ffi.string(resource[0]).decode("UTF-8") 78 79 # Preferred MemType 80 def get_preferred_memtype(self): 81 """Return Ceed preferred memory type. 82 83 Returns: 84 memtype: Ceed preferred memory type""" 85 86 # libCEED call 87 memtype = ffi.new("CeedMemType *", MEM_HOST) 88 lib.CeedGetPreferredMemType(self._pointer[0], memtype) 89 90 return memtype[0] 91 92 # CeedVector 93 def Vector(self, size): 94 """Ceed Vector: storing and manipulating vectors. 95 96 Args: 97 size: length of vector 98 99 Returns: 100 vector: Ceed Vector""" 101 102 return Vector(self, size) 103 104 # CeedElemRestriction 105 def ElemRestriction(self, nelem, elemsize, nnodes, ncomp, indices, 106 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES, 107 imode=lib.CEED_NONINTERLACED): 108 """Ceed ElemRestriction: restriction from local vectors to elements. 109 110 Args: 111 nelem: number of elements described in the indices array 112 elemsize: size (number of nodes) per element 113 nnodes: the number of nodes in the local vector. The input Ceed Vector 114 to which the restriction will be applied is of size 115 (nnodes * ncomp). This size may include data 116 used by other Ceed ElemRestriction objects describing 117 different types of elements. 118 ncomp: number of field components per interpolation node (1 for scalar fields) 119 *indices: Numpy array of shape [nelem, elemsize]. Row i holds the 120 ordered list of the indices (into the input Ceed Vector) 121 for the unknowns corresponding to element i, where 122 0 <= i < nelem. All indices must be in the range 123 [0, nnodes - 1]. 124 **memtype: memory type of the indices array, default CEED_MEM_HOST 125 **cmode: copy mode for the indices array, default CEED_COPY_VALUES 126 **imode: ordering of the ncomp components, i.e. it specifies 127 the ordering of the components of the local vector used 128 by this CeedElemRestriction. CEED_NONINTERLACED indicates 129 the component is the outermost index and CEED_INTERLACED 130 indicates the component is the innermost index in 131 ordering of the local vector, default CEED_NONINTERLACED 132 133 Returns: 134 elemrestriction: Ceed ElemRestiction""" 135 136 return ElemRestriction(self, nelem, elemsize, nnodes, ncomp, indices, 137 memtype=memtype, cmode=cmode, imode=imode) 138 139 def StridedElemRestriction(self, nelem, elemsize, nnodes, ncomp, strides): 140 """Ceed Identity ElemRestriction: strided restriction from local vectors to elements. 141 142 Args: 143 nelem: number of elements described in the indices array 144 elemsize: size (number of nodes) per element 145 nnodes: the number of nodes in the local vector. The input Ceed Vector 146 to which the restriction will be applied is of size 147 (nnodes * ncomp). This size may include data 148 used by other Ceed ElemRestriction objects describing 149 different types of elements. 150 ncomp: number of field components per interpolation node 151 (1 for scalar fields) 152 *strides: Array for strides between [nodes, components, elements]. 153 The data for node i, component j, element k in the 154 L-vector is given by 155 i*strides[0] + j*strides[1] + k*strides[2] 156 157 Returns: 158 elemrestriction: Ceed Strided ElemRestiction""" 159 160 return StridedElemRestriction( 161 self, nelem, elemsize, nnodes, ncomp, strides) 162 163 def BlockedElemRestriction(self, nelem, elemsize, blksize, nnodes, ncomp, 164 indices, memtype=lib.CEED_MEM_HOST, 165 cmode=lib.CEED_COPY_VALUES, 166 imode=lib.CEED_NONINTERLACED): 167 """Ceed Blocked ElemRestriction: blocked restriction from local vectors to elements. 168 169 Args: 170 nelem: number of elements described in the indices array 171 elemsize: size (number of nodes) per element 172 blksize: number of elements in a block 173 nnodes: the number of nodes in the local vector. The input Ceed Vector 174 to which the restriction will be applied is of size 175 (nnodes * ncomp). This size may include data 176 used by other Ceed ElemRestriction objects describing 177 different types of elements. 178 ncomp: number of field components per interpolation node (1 for scalar fields) 179 *indices: Numpy array of shape [nelem, elemsize]. Row i holds the 180 ordered list of the indices (into the input Ceed Vector) 181 for the unknowns corresponding to element i, where 182 0 <= i < nelem. All indices must be in the range 183 [0, nnodes - 1]. The backend will permute and pad this 184 array to the desired ordering for the blocksize, which is 185 typically given by the backend. The default reordering is 186 to interlace elements. 187 **memtype: memory type of the indices array, default CEED_MEM_HOST 188 **cmode: copy mode for the indices array, default CEED_COPY_VALUES 189 **imode: ordering of the ncomp components, i.e. it specifies 190 the ordering of the components of the local vector used 191 by this CeedElemRestriction. CEED_NONINTERLACED indicates 192 the component is the outermost index and CEED_INTERLACED 193 indicates the component is the innermost index in 194 ordering of the local vector, default CEED_NONINTERLACED 195 196 Returns: 197 elemrestriction: Ceed Blocked ElemRestiction""" 198 199 return BlockedElemRestriction(self, nelem, elemsize, blksize, nnodes, 200 ncomp, indices, memtype=memtype, 201 cmode=cmode, imode=imode) 202 203 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, nnodes, 204 ncomp, strides): 205 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction from local vectors to elements. 206 207 Args: 208 nelem: number of elements described in the indices array 209 elemsize: size (number of nodes) per element 210 blksize: number of elements in a block 211 nnodes: the number of nodes in the local vector. The input Ceed Vector 212 to which the restriction will be applied is of size 213 (nnodes * ncomp). This size may include data 214 used by other Ceed ElemRestriction objects describing 215 different types of elements. 216 ncomp: number of field components per interpolation node 217 (1 for scalar fields) 218 *strides: Array for strides between [nodes, components, elements]. 219 The data for node i, component j, element k in the 220 L-vector is given by 221 i*strides[0] + j*strides[1] + k*strides[2] 222 223 Returns: 224 elemrestriction: Ceed Strided ElemRestiction""" 225 226 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 227 nnodes, ncomp, strides) 228 229 # CeedBasis 230 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 231 qref1d, qweight1d): 232 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 233 H^1 discretizations. 234 235 Args: 236 dim: topological dimension 237 ncomp: number of field components (1 for scalar fields) 238 P1d: number of nodes in one dimension 239 Q1d: number of quadrature points in one dimension 240 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix expressing the 241 values of nodal basis functions at quadrature points 242 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix expressing the 243 derivatives of nodal basis functions at quadrature points 244 *qref1d: Numpy array of length Q1d holding the locations of quadrature points 245 on the 1D reference element [-1, 1] 246 *qweight1d: Numpy array of length Q1d holding the quadrature weights on the 247 reference element 248 249 Returns: 250 basis: Ceed Basis""" 251 252 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 253 qref1d, qweight1d) 254 255 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 256 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange basis 257 objects for H^1 discretizations. 258 259 Args: 260 dim: topological dimension 261 ncomp: number of field components (1 for scalar fields) 262 P: number of Gauss-Lobatto nodes in one dimension. The 263 polynomial degree of the resulting Q_k element is k=P-1. 264 Q: number of quadrature points in one dimension 265 qmode: distribution of the Q quadrature points (affects order of 266 accuracy for the quadrature) 267 268 Returns: 269 basis: Ceed Basis""" 270 271 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 272 273 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 274 """Ceed H1 Basis: finite element non tensor-product basis for H^1 discretizations. 275 276 Args: 277 topo: topology of the element, e.g. hypercube, simplex, etc 278 ncomp: number of field components (1 for scalar fields) 279 nnodes: total number of nodes 280 nqpts: total number of quadrature points 281 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix expressing 282 the values of nodal basis functions at quadrature points 283 *grad: Numpy array holding the row-major (nqpts * dim * nnodes) matrix 284 expressing the derivatives of nodal basis functions at 285 quadrature points 286 *qref: Numpy array of length (nqpts * dim) holding the locations of quadrature 287 points on the reference element [-1, 1] 288 *qweight: Numpy array of length nnodes holding the quadrature weights on the 289 reference element 290 291 Returns: 292 basis: Ceed Basis""" 293 294 return BasisH1(self, topo, ncomp, nnodes, nqpts, 295 interp, grad, qref, qweight) 296 297 # CeedQFunction 298 def QFunction(self, vlength, f, source): 299 """Ceed QFunction: point-wise operation at quadrature points for evaluating 300 volumetric terms. 301 302 Args: 303 vlength: vector length. Caller must ensure that number of quadrature 304 points is a multiple of vlength 305 f: ctypes function pointer to evaluate action at quadrature points 306 source: absolute path to source of QFunction, 307 "\\abs_path\\file.h:function_name 308 309 Returns: 310 qfunction: Ceed QFunction""" 311 312 return QFunction(self, vlength, f, source) 313 314 def QFunctionByName(self, name): 315 """Ceed QFunction By Name: point-wise operation at quadrature points 316 from a given gallery, for evaluating volumetric terms. 317 318 Args: 319 name: name of QFunction to use from gallery 320 321 Returns: 322 qfunction: Ceed QFunction By Name""" 323 324 return QFunctionByName(self, name) 325 326 def IdentityQFunction(self, size, inmode, outmode): 327 """Ceed Idenity QFunction: identity qfunction operation. 328 329 Args: 330 size: size of the qfunction fields 331 **inmode: CeedEvalMode for input to Ceed QFunction 332 **outmode: CeedEvalMode for output to Ceed QFunction 333 334 Returns: 335 qfunction: Ceed Identity QFunction""" 336 337 return IdentityQFunction(self, size, inmode, outmode) 338 339 # CeedOperator 340 def Operator(self, qf, dqf=None, qdfT=None): 341 """Ceed Operator: composed FE-type operations on vectors. 342 343 Args: 344 qf: QFunction defining the action of the operator at quadrature points 345 **dqf: QFunction defining the action of the Jacobian, default None 346 **dqfT: QFunction defining the action of the transpose of the Jacobian, 347 default None 348 349 Returns: 350 operator: Ceed Operator""" 351 352 return Operator(self, qf, dqf, qdfT) 353 354 def CompositeOperator(self): 355 """Composite Ceed Operator: composition of multiple CeedOperators. 356 357 Returns: 358 operator: Ceed Composite Operator""" 359 360 return CompositeOperator(self) 361 362 # Destructor 363 def __del__(self): 364 # libCEED call 365 lib.CeedDestroy(self._pointer) 366 367# ------------------------------------------------------------------------------ 368