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 tempfile 19import numpy as np 20from abc import ABC 21from .ceed_constants import TRANSPOSE, NOTRANSPOSE 22 23# ------------------------------------------------------------------------------ 24 25 26class Basis(ABC): 27 """Ceed Basis: finite element basis objects.""" 28 29 # Representation 30 def __repr__(self): 31 return "<CeedBasis instance at " + hex(id(self)) + ">" 32 33 # String conversion for print() to stdout 34 def __str__(self): 35 """View a Basis via print().""" 36 37 # libCEED call 38 with tempfile.NamedTemporaryFile() as key_file: 39 with open(key_file.name, 'r+') as stream_file: 40 stream = ffi.cast("FILE *", stream_file) 41 42 err_code = lib.CeedBasisView(self._pointer[0], stream) 43 self._ceed._check_error(err_code) 44 45 stream_file.seek(0) 46 out_string = stream_file.read() 47 48 return out_string 49 50 # Apply Basis 51 def apply(self, nelem, emode, u, v, tmode=NOTRANSPOSE): 52 """Apply basis evaluation from nodes to quadrature points or vice versa. 53 54 Args: 55 nelem: the number of elements to apply the basis evaluation to; 56 the backend will specify the ordering in a 57 BlockedElemRestriction 58 emode: basis evaluation mode 59 u: input vector 60 v: output vector 61 **tmode: CEED_NOTRANSPOSE to evaluate from nodes to quadrature 62 points, CEED_TRANSPOSE to apply the transpose, mapping 63 from quadrature points to nodes; default CEED_NOTRANSPOSE""" 64 65 # libCEED call 66 err_code = lib.CeedBasisApply(self._pointer[0], nelem, tmode, emode, 67 u._pointer[0], v._pointer[0]) 68 self._ceed._check_error(err_code) 69 70 # Transpose a Basis 71 @property 72 def T(self): 73 """Transpose a Basis.""" 74 75 return TransposeBasis(self) 76 77 # Transpose a Basis 78 @property 79 def transpose(self): 80 """Transpose a Basis.""" 81 82 return TransposeBasis(self) 83 84 # Get number of nodes 85 def get_num_nodes(self): 86 """Get total number of nodes (in dim dimensions) of a Basis. 87 88 Returns: 89 num_nodes: total number of nodes""" 90 91 # Setup argument 92 p_pointer = ffi.new("CeedInt *") 93 94 # libCEED call 95 err_code = lib.CeedBasisGetNumNodes(self._pointer[0], p_pointer) 96 self._ceed._check_error(err_code) 97 98 return p_pointer[0] 99 100 # Get number of quadrature points 101 def get_num_quadrature_points(self): 102 """Get total number of quadrature points (in dim dimensions) of a Basis. 103 104 Returns: 105 num_qpts: total number of quadrature points""" 106 107 # Setup argument 108 q_pointer = ffi.new("CeedInt *") 109 110 # libCEED call 111 err_code = lib.CeedBasisGetNumQuadraturePoints( 112 self._pointer[0], q_pointer) 113 self._ceed._check_error(err_code) 114 115 return q_pointer[0] 116 117 # Destructor 118 def __del__(self): 119 # libCEED call 120 err_code = lib.CeedBasisDestroy(self._pointer) 121 self._ceed._check_error(err_code) 122 123# ------------------------------------------------------------------------------ 124 125 126class BasisTensorH1(Basis): 127 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 128 H^1 discretizations.""" 129 130 # Constructor 131 def __init__(self, ceed, dim, ncomp, P1d, Q1d, interp1d, grad1d, 132 qref1d, qweight1d): 133 134 # Setup arguments 135 self._pointer = ffi.new("CeedBasis *") 136 137 self._ceed = ceed 138 139 interp1d_pointer = ffi.new("CeedScalar *") 140 interp1d_pointer = ffi.cast( 141 "CeedScalar *", 142 interp1d.__array_interface__['data'][0]) 143 144 grad1d_pointer = ffi.new("CeedScalar *") 145 grad1d_pointer = ffi.cast( 146 "CeedScalar *", 147 grad1d.__array_interface__['data'][0]) 148 149 qref1d_pointer = ffi.new("CeedScalar *") 150 qref1d_pointer = ffi.cast( 151 "CeedScalar *", 152 qref1d.__array_interface__['data'][0]) 153 154 qweight1d_pointer = ffi.new("CeedScalar *") 155 qweight1d_pointer = ffi.cast( 156 "CeedScalar *", 157 qweight1d.__array_interface__['data'][0]) 158 159 # libCEED call 160 err_code = lib.CeedBasisCreateTensorH1(self._ceed._pointer[0], dim, ncomp, 161 P1d, Q1d, interp1d_pointer, 162 grad1d_pointer, qref1d_pointer, 163 qweight1d_pointer, self._pointer) 164 self._ceed._check_error(err_code) 165 166 # 167 def get_interp1d(self): 168 """Return 1D interpolation matrix of a tensor product Basis. 169 170 Returns: 171 *array: Numpy array""" 172 173 # Compute the length of the array 174 nnodes_pointer = ffi.new("CeedInt *") 175 lib.CeedBasisGetNumNodes1D(self._pointer[0], nnodes_pointer) 176 nqpts_pointer = ffi.new("CeedInt *") 177 lib.CeedBasisGetNumQuadraturePoints1D(self._pointer[0], nqpts_pointer) 178 length = nnodes_pointer[0] * nqpts_pointer[0] 179 180 # Setup the pointer's pointer 181 array_pointer = ffi.new("CeedScalar **") 182 183 # libCEED call 184 lib.CeedBasisGetInterp1D(self._pointer[0], array_pointer) 185 186 # Return array created from buffer 187 # Create buffer object from returned pointer 188 buff = ffi.buffer( 189 array_pointer[0], 190 ffi.sizeof("CeedScalar") * 191 length) 192 # return read only Numpy array 193 ret = np.frombuffer(buff, dtype="float64") 194 ret.flags['WRITEABLE'] = False 195 return ret 196 197 198# ------------------------------------------------------------------------------ 199 200 201class BasisTensorH1Lagrange(BasisTensorH1): 202 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange basis 203 objects for H^1 discretizations.""" 204 205 # Constructor 206 def __init__(self, ceed, dim, ncomp, P, Q, qmode): 207 208 # Setup arguments 209 self._pointer = ffi.new("CeedBasis *") 210 211 self._ceed = ceed 212 213 # libCEED call 214 err_code = lib.CeedBasisCreateTensorH1Lagrange(self._ceed._pointer[0], dim, 215 ncomp, P, Q, qmode, self._pointer) 216 self._ceed._check_error(err_code) 217 218# ------------------------------------------------------------------------------ 219 220 221class BasisH1(Basis): 222 """Ceed H1 Basis: finite element non tensor-product basis for H^1 discretizations.""" 223 224 # Constructor 225 def __init__(self, ceed, topo, ncomp, nnodes, 226 nqpts, interp, grad, qref, qweight): 227 228 # Setup arguments 229 self._pointer = ffi.new("CeedBasis *") 230 231 self._ceed = ceed 232 233 interp_pointer = ffi.new("CeedScalar *") 234 interp_pointer = ffi.cast( 235 "CeedScalar *", 236 interp.__array_interface__['data'][0]) 237 238 grad_pointer = ffi.new("CeedScalar *") 239 grad_pointer = ffi.cast( 240 "CeedScalar *", 241 grad.__array_interface__['data'][0]) 242 243 qref_pointer = ffi.new("CeedScalar *") 244 qref_pointer = ffi.cast( 245 "CeedScalar *", 246 qref.__array_interface__['data'][0]) 247 248 qweight_pointer = ffi.new("CeedScalar *") 249 qweight_pointer = ffi.cast( 250 "CeedScalar *", 251 qweight.__array_interface__['data'][0]) 252 253 # libCEED call 254 err_code = lib.CeedBasisCreateH1(self._ceed._pointer[0], topo, ncomp, 255 nnodes, nqpts, interp_pointer, 256 grad_pointer, qref_pointer, 257 qweight_pointer, self._pointer) 258 259# ------------------------------------------------------------------------------ 260 261 262class TransposeBasis(): 263 """Transpose Ceed Basis: transpose of finite element tensor-product basis objects.""" 264 265 # Attributes 266 _basis = None 267 268 # Constructor 269 def __init__(self, basis): 270 271 # Reference basis 272 self._basis = basis 273 274 # Representation 275 def __repr__(self): 276 return "<Transpose CeedBasis instance at " + hex(id(self)) + ">" 277 278 # Apply Transpose Basis 279 def apply(self, nelem, emode, u, v): 280 """Apply basis evaluation from quadrature points to nodes. 281 282 Args: 283 nelem: the number of elements to apply the basis evaluation to; 284 the backend will specify the ordering in a 285 Blocked ElemRestriction 286 **emode: basis evaluation mode 287 u: input vector 288 v: output vector""" 289 290 # libCEED call 291 self._basis.apply(nelem, emode, u, v, tmode=TRANSPOSE) 292 293# ------------------------------------------------------------------------------ 294