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