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, ncomp, compstride, lsize, offsets, 106 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES): 107 """Ceed ElemRestriction: restriction from local vectors to elements. 108 109 Args: 110 nelem: number of elements described by the restriction 111 elemsize: size (number of nodes) per element 112 ncomp: number of field components per interpolation node 113 (1 for scalar fields) 114 compstride: Stride between components for the same L-vector "node". 115 Data for node i, component k can be found in the 116 L-vector at index [offsets[i] + k*compstride]. 117 lsize: The size of the L-vector. This vector may be larger than 118 the elements and fields given by this restriction. 119 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 120 ordered list of the offsets (into the input Ceed Vector) 121 for the unknowns corresponding to element i, where 122 0 <= i < nelem. All offsets must be in the range 123 [0, lsize - 1]. 124 **memtype: memory type of the offsets array, default CEED_MEM_HOST 125 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 126 127 Returns: 128 elemrestriction: Ceed ElemRestiction""" 129 130 return ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, 131 offsets, memtype=memtype, cmode=cmode) 132 133 def StridedElemRestriction(self, nelem, elemsize, ncomp, lsize, strides): 134 """Ceed Identity ElemRestriction: strided restriction from local vectors 135 to elements. 136 137 Args: 138 nelem: number of elements described by the restriction 139 elemsize: size (number of nodes) per element 140 ncomp: number of field components per interpolation node 141 (1 for scalar fields) 142 lsize: The size of the L-vector. This vector may be larger than 143 the elements and fields given by this restriction. 144 *strides: Array for strides between [nodes, components, elements]. 145 The data for node i, component j, element k in the 146 L-vector is given by 147 i*strides[0] + j*strides[1] + k*strides[2] 148 149 Returns: 150 elemrestriction: Ceed Strided ElemRestiction""" 151 152 return StridedElemRestriction( 153 self, nelem, elemsize, ncomp, lsize, strides) 154 155 def BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, compstride, 156 lsize, offsets, memtype=lib.CEED_MEM_HOST, 157 cmode=lib.CEED_COPY_VALUES): 158 """Ceed Blocked ElemRestriction: blocked restriction from local vectors 159 to elements. 160 161 Args: 162 nelem: number of elements described by the restriction 163 elemsize: size (number of nodes) per element 164 blksize: number of elements in a block 165 ncomp: number of field components per interpolation node 166 (1 for scalar fields) 167 lsize: The size of the L-vector. This vector may be larger than 168 the elements and fields given by this restriction. 169 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 170 ordered list of the offsets (into the input Ceed Vector) 171 for the unknowns corresponding to element i, where 172 0 <= i < nelem. All offsets must be in the range 173 [0, lsize - 1]. The backend will permute and pad this 174 array to the desired ordering for the blocksize, which is 175 typically given by the backend. The default reordering is 176 to interlace elements. 177 **memtype: memory type of the offsets array, default CEED_MEM_HOST 178 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 179 180 Returns: 181 elemrestriction: Ceed Blocked ElemRestiction""" 182 183 return BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, 184 compstride, lsize, offsets, 185 memtype=memtype, cmode=cmode) 186 187 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, ncomp, 188 lsize, strides): 189 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction 190 from local vectors to elements. 191 192 Args: 193 nelem: number of elements described in the restriction 194 elemsize: size (number of nodes) per element 195 blksize: number of elements in a block 196 ncomp: number of field components per interpolation node 197 (1 for scalar fields) 198 lsize: The size of the L-vector. This vector may be larger than 199 the elements and fields given by this restriction. 200 *strides: Array for strides between [nodes, components, elements]. 201 The data for node i, component j, element k in the 202 L-vector is given by 203 i*strides[0] + j*strides[1] + k*strides[2] 204 205 Returns: 206 elemrestriction: Ceed Strided ElemRestiction""" 207 208 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 209 ncomp, lsize, strides) 210 211 # CeedBasis 212 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 213 qref1d, qweight1d): 214 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 215 H^1 discretizations. 216 217 Args: 218 dim: topological dimension 219 ncomp: number of field components (1 for scalar fields) 220 P1d: number of nodes in one dimension 221 Q1d: number of quadrature points in one dimension 222 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix 223 expressing the values of nodal basis functions at 224 quadrature points 225 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix 226 expressing the derivatives of nodal basis functions at 227 quadrature points 228 *qref1d: Numpy array of length Q1d holding the locations of 229 quadrature points on the 1D reference element [-1, 1] 230 *qweight1d: Numpy array of length Q1d holding the quadrature 231 weights on the reference element 232 233 Returns: 234 basis: Ceed Basis""" 235 236 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 237 qref1d, qweight1d) 238 239 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 240 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange 241 basis objects for H^1 discretizations. 242 243 Args: 244 dim: topological dimension 245 ncomp: number of field components (1 for scalar fields) 246 P: number of Gauss-Lobatto nodes in one dimension. The 247 polynomial degree of the resulting Q_k element is k=P-1. 248 Q: number of quadrature points in one dimension 249 qmode: distribution of the Q quadrature points (affects order of 250 accuracy for the quadrature) 251 252 Returns: 253 basis: Ceed Basis""" 254 255 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 256 257 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 258 """Ceed H1 Basis: finite element non tensor-product basis for H^1 259 discretizations. 260 261 Args: 262 topo: topology of the element, e.g. hypercube, simplex, etc 263 ncomp: number of field components (1 for scalar fields) 264 nnodes: total number of nodes 265 nqpts: total number of quadrature points 266 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix 267 expressing the values of nodal basis functions at 268 quadrature points 269 *grad: Numpy array holding the row-major (nqpts * dim * nnodes) 270 matrix expressing the derivatives of nodal basis functions 271 at quadrature points 272 *qref: Numpy array of length (nqpts * dim) holding the locations of 273 quadrature points on the reference element [-1, 1] 274 *qweight: Numpy array of length nnodes holding the quadrature 275 weights on the reference element 276 277 Returns: 278 basis: Ceed Basis""" 279 280 return BasisH1(self, topo, ncomp, nnodes, nqpts, 281 interp, grad, qref, qweight) 282 283 # CeedQFunction 284 def QFunction(self, vlength, f, source): 285 """Ceed QFunction: point-wise operation at quadrature points for 286 evaluating volumetric terms. 287 288 Args: 289 vlength: vector length. Caller must ensure that number of quadrature 290 points is a multiple of vlength 291 f: ctypes function pointer to evaluate action at quadrature points 292 source: absolute path to source of QFunction, 293 "\\abs_path\\file.h:function_name 294 295 Returns: 296 qfunction: Ceed QFunction""" 297 298 return QFunction(self, vlength, f, source) 299 300 def QFunctionByName(self, name): 301 """Ceed QFunction By Name: point-wise operation at quadrature points 302 from a given gallery, for evaluating volumetric terms. 303 304 Args: 305 name: name of QFunction to use from gallery 306 307 Returns: 308 qfunction: Ceed QFunction By Name""" 309 310 return QFunctionByName(self, name) 311 312 def IdentityQFunction(self, size, inmode, outmode): 313 """Ceed Idenity QFunction: identity qfunction operation. 314 315 Args: 316 size: size of the qfunction fields 317 **inmode: CeedEvalMode for input to Ceed QFunction 318 **outmode: CeedEvalMode for output to Ceed QFunction 319 320 Returns: 321 qfunction: Ceed Identity QFunction""" 322 323 return IdentityQFunction(self, size, inmode, outmode) 324 325 # CeedOperator 326 def Operator(self, qf, dqf=None, qdfT=None): 327 """Ceed Operator: composed FE-type operations on vectors. 328 329 Args: 330 qf: QFunction defining the action of the operator at quadrature 331 points 332 **dqf: QFunction defining the action of the Jacobian, default None 333 **dqfT: QFunction defining the action of the transpose of the 334 Jacobian, default None 335 336 Returns: 337 operator: Ceed Operator""" 338 339 return Operator(self, qf, dqf, qdfT) 340 341 def CompositeOperator(self): 342 """Composite Ceed Operator: composition of multiple CeedOperators. 343 344 Returns: 345 operator: Ceed Composite Operator""" 346 347 return CompositeOperator(self) 348 349 # Destructor 350 def __del__(self): 351 # libCEED call 352 lib.CeedDestroy(self._pointer) 353 354# ------------------------------------------------------------------------------ 355