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 17 //! A Ceed Basis defines the discrete finite element basis and associated 18 //! quadrature rule. 19 20 use crate::prelude::*; 21 22 // ----------------------------------------------------------------------------- 23 // CeedBasis option 24 // ----------------------------------------------------------------------------- 25 #[derive(Clone, Copy)] 26 pub enum BasisOpt<'a> { 27 Some(&'a Basis<'a>), 28 Collocated, 29 } 30 /// Construct a BasisOpt reference from a Basis reference 31 impl<'a> From<&'a Basis<'_>> for BasisOpt<'a> { 32 fn from(basis: &'a Basis) -> Self { 33 debug_assert!(basis.ptr != unsafe { bind_ceed::CEED_BASIS_COLLOCATED }); 34 Self::Some(basis) 35 } 36 } 37 impl<'a> BasisOpt<'a> { 38 /// Transform a Rust libCEED BasisOpt into C libCEED CeedBasis 39 pub(crate) fn to_raw(self) -> bind_ceed::CeedBasis { 40 match self { 41 Self::Some(basis) => basis.ptr, 42 Self::Collocated => unsafe { bind_ceed::CEED_BASIS_COLLOCATED }, 43 } 44 } 45 } 46 47 // ----------------------------------------------------------------------------- 48 // CeedBasis context wrapper 49 // ----------------------------------------------------------------------------- 50 pub struct Basis<'a> { 51 ceed: &'a crate::Ceed, 52 pub(crate) ptr: bind_ceed::CeedBasis, 53 } 54 55 // ----------------------------------------------------------------------------- 56 // Destructor 57 // ----------------------------------------------------------------------------- 58 impl<'a> Drop for Basis<'a> { 59 fn drop(&mut self) { 60 unsafe { 61 if self.ptr != bind_ceed::CEED_BASIS_COLLOCATED { 62 bind_ceed::CeedBasisDestroy(&mut self.ptr); 63 } 64 } 65 } 66 } 67 68 // ----------------------------------------------------------------------------- 69 // Display 70 // ----------------------------------------------------------------------------- 71 impl<'a> fmt::Display for Basis<'a> { 72 /// View a Basis 73 /// 74 /// ``` 75 /// # use libceed::prelude::*; 76 /// # let ceed = libceed::Ceed::default_init(); 77 /// let b = ceed 78 /// .basis_tensor_H1_Lagrange(1, 2, 3, 4, QuadMode::Gauss) 79 /// .unwrap(); 80 /// println!("{}", b); 81 /// ``` 82 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 83 let mut ptr = std::ptr::null_mut(); 84 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 85 let cstring = unsafe { 86 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 87 bind_ceed::CeedBasisView(self.ptr, file); 88 bind_ceed::fclose(file); 89 CString::from_raw(ptr) 90 }; 91 cstring.to_string_lossy().fmt(f) 92 } 93 } 94 95 // ----------------------------------------------------------------------------- 96 // Implementations 97 // ----------------------------------------------------------------------------- 98 impl<'a> Basis<'a> { 99 // Constructors 100 pub fn create_tensor_H1( 101 ceed: &'a crate::Ceed, 102 dim: usize, 103 ncomp: usize, 104 P1d: usize, 105 Q1d: usize, 106 interp1d: &[crate::Scalar], 107 grad1d: &[crate::Scalar], 108 qref1d: &[crate::Scalar], 109 qweight1d: &[crate::Scalar], 110 ) -> crate::Result<Self> { 111 let mut ptr = std::ptr::null_mut(); 112 let (dim, ncomp, P1d, Q1d) = ( 113 i32::try_from(dim).unwrap(), 114 i32::try_from(ncomp).unwrap(), 115 i32::try_from(P1d).unwrap(), 116 i32::try_from(Q1d).unwrap(), 117 ); 118 let ierr = unsafe { 119 bind_ceed::CeedBasisCreateTensorH1( 120 ceed.ptr, 121 dim, 122 ncomp, 123 P1d, 124 Q1d, 125 interp1d.as_ptr(), 126 grad1d.as_ptr(), 127 qref1d.as_ptr(), 128 qweight1d.as_ptr(), 129 &mut ptr, 130 ) 131 }; 132 ceed.check_error(ierr)?; 133 Ok(Self { ceed, ptr }) 134 } 135 136 pub fn create_tensor_H1_Lagrange( 137 ceed: &'a crate::Ceed, 138 dim: usize, 139 ncomp: usize, 140 P: usize, 141 Q: usize, 142 qmode: crate::QuadMode, 143 ) -> crate::Result<Self> { 144 let mut ptr = std::ptr::null_mut(); 145 let (dim, ncomp, P, Q, qmode) = ( 146 i32::try_from(dim).unwrap(), 147 i32::try_from(ncomp).unwrap(), 148 i32::try_from(P).unwrap(), 149 i32::try_from(Q).unwrap(), 150 qmode as bind_ceed::CeedQuadMode, 151 ); 152 let ierr = unsafe { 153 bind_ceed::CeedBasisCreateTensorH1Lagrange(ceed.ptr, dim, ncomp, P, Q, qmode, &mut ptr) 154 }; 155 ceed.check_error(ierr)?; 156 Ok(Self { ceed, ptr }) 157 } 158 159 pub fn create_H1( 160 ceed: &'a crate::Ceed, 161 topo: crate::ElemTopology, 162 ncomp: usize, 163 nnodes: usize, 164 nqpts: usize, 165 interp: &[crate::Scalar], 166 grad: &[crate::Scalar], 167 qref: &[crate::Scalar], 168 qweight: &[crate::Scalar], 169 ) -> crate::Result<Self> { 170 let mut ptr = std::ptr::null_mut(); 171 let (topo, ncomp, nnodes, nqpts) = ( 172 topo as bind_ceed::CeedElemTopology, 173 i32::try_from(ncomp).unwrap(), 174 i32::try_from(nnodes).unwrap(), 175 i32::try_from(nqpts).unwrap(), 176 ); 177 let ierr = unsafe { 178 bind_ceed::CeedBasisCreateH1( 179 ceed.ptr, 180 topo, 181 ncomp, 182 nnodes, 183 nqpts, 184 interp.as_ptr(), 185 grad.as_ptr(), 186 qref.as_ptr(), 187 qweight.as_ptr(), 188 &mut ptr, 189 ) 190 }; 191 ceed.check_error(ierr)?; 192 Ok(Self { ceed, ptr }) 193 } 194 195 /// Apply basis evaluation from nodes to quadrature points or vice versa 196 /// 197 /// * `nelem` - The number of elements to apply the basis evaluation to 198 /// * `tmode` - `TrasposeMode::NoTranspose` to evaluate from nodes to 199 /// quadrature points, `TransposeMode::Transpose` to apply the 200 /// transpose, mapping from quadrature points to nodes 201 /// * `emode` - `EvalMode::None` to use values directly, `EvalMode::Interp` 202 /// to use interpolated values, `EvalMode::Grad` to use 203 /// gradients, `EvalMode::Weight` to use quadrature weights 204 /// * `u` - Input Vector 205 /// * `v` - Output Vector 206 /// 207 /// ``` 208 /// # use libceed::prelude::*; 209 /// # let ceed = libceed::Ceed::default_init(); 210 /// const Q: usize = 6; 211 /// let bu = ceed 212 /// .basis_tensor_H1_Lagrange(1, 1, Q, Q, QuadMode::GaussLobatto) 213 /// .unwrap(); 214 /// let bx = ceed 215 /// .basis_tensor_H1_Lagrange(1, 1, 2, Q, QuadMode::Gauss) 216 /// .unwrap(); 217 /// 218 /// let x_corners = ceed.vector_from_slice(&[-1., 1.]).unwrap(); 219 /// let mut x_qpts = ceed.vector(Q).unwrap(); 220 /// let mut x_nodes = ceed.vector(Q).unwrap(); 221 /// bx.apply( 222 /// 1, 223 /// TransposeMode::NoTranspose, 224 /// EvalMode::Interp, 225 /// &x_corners, 226 /// &mut x_nodes, 227 /// ); 228 /// bu.apply( 229 /// 1, 230 /// TransposeMode::NoTranspose, 231 /// EvalMode::Interp, 232 /// &x_nodes, 233 /// &mut x_qpts, 234 /// ); 235 /// 236 /// // Create function x^3 + 1 on Gauss Lobatto points 237 /// let mut u_arr = [0.; Q]; 238 /// u_arr 239 /// .iter_mut() 240 /// .zip(x_nodes.view().iter()) 241 /// .for_each(|(u, x)| *u = x * x * x + 1.); 242 /// let u = ceed.vector_from_slice(&u_arr).unwrap(); 243 /// 244 /// // Map function to Gauss points 245 /// let mut v = ceed.vector(Q).unwrap(); 246 /// v.set_value(0.); 247 /// bu.apply(1, TransposeMode::NoTranspose, EvalMode::Interp, &u, &mut v) 248 /// .unwrap(); 249 /// 250 /// // Verify results 251 /// v.view() 252 /// .iter() 253 /// .zip(x_qpts.view().iter()) 254 /// .for_each(|(v, x)| { 255 /// let true_value = x * x * x + 1.; 256 /// assert!( 257 /// (*v - true_value).abs() < 10.0 * libceed::EPSILON, 258 /// "Incorrect basis application" 259 /// ); 260 /// }); 261 /// ``` 262 pub fn apply( 263 &self, 264 nelem: usize, 265 tmode: TransposeMode, 266 emode: EvalMode, 267 u: &Vector, 268 v: &mut Vector, 269 ) -> crate::Result<i32> { 270 let (nelem, tmode, emode) = ( 271 i32::try_from(nelem).unwrap(), 272 tmode as bind_ceed::CeedTransposeMode, 273 emode as bind_ceed::CeedEvalMode, 274 ); 275 let ierr = 276 unsafe { bind_ceed::CeedBasisApply(self.ptr, nelem, tmode, emode, u.ptr, v.ptr) }; 277 self.ceed.check_error(ierr) 278 } 279 280 /// Returns the dimension for given CeedBasis 281 /// 282 /// ``` 283 /// # use libceed::prelude::*; 284 /// # let ceed = libceed::Ceed::default_init(); 285 /// let dim = 2; 286 /// let b = ceed 287 /// .basis_tensor_H1_Lagrange(dim, 1, 3, 4, QuadMode::Gauss) 288 /// .unwrap(); 289 /// 290 /// let d = b.dimension(); 291 /// assert_eq!(d, dim, "Incorrect dimension"); 292 /// ``` 293 pub fn dimension(&self) -> usize { 294 let mut dim = 0; 295 unsafe { bind_ceed::CeedBasisGetDimension(self.ptr, &mut dim) }; 296 usize::try_from(dim).unwrap() 297 } 298 299 /// Returns number of components for given CeedBasis 300 /// 301 /// ``` 302 /// # use libceed::prelude::*; 303 /// # let ceed = libceed::Ceed::default_init(); 304 /// let ncomp = 2; 305 /// let b = ceed 306 /// .basis_tensor_H1_Lagrange(1, ncomp, 3, 4, QuadMode::Gauss) 307 /// .unwrap(); 308 /// 309 /// let n = b.num_components(); 310 /// assert_eq!(n, ncomp, "Incorrect number of components"); 311 /// ``` 312 pub fn num_components(&self) -> usize { 313 let mut ncomp = 0; 314 unsafe { bind_ceed::CeedBasisGetNumComponents(self.ptr, &mut ncomp) }; 315 usize::try_from(ncomp).unwrap() 316 } 317 318 /// Returns total number of nodes (in dim dimensions) of a CeedBasis 319 /// 320 /// ``` 321 /// # use libceed::prelude::*; 322 /// # let ceed = libceed::Ceed::default_init(); 323 /// let p = 3; 324 /// let b = ceed 325 /// .basis_tensor_H1_Lagrange(2, 1, p, 4, QuadMode::Gauss) 326 /// .unwrap(); 327 /// 328 /// let nnodes = b.num_nodes(); 329 /// assert_eq!(nnodes, p * p, "Incorrect number of nodes"); 330 /// ``` 331 pub fn num_nodes(&self) -> usize { 332 let mut nnodes = 0; 333 unsafe { bind_ceed::CeedBasisGetNumNodes(self.ptr, &mut nnodes) }; 334 usize::try_from(nnodes).unwrap() 335 } 336 337 /// Returns total number of quadrature points (in dim dimensions) of a 338 /// CeedBasis 339 /// 340 /// ``` 341 /// # use libceed::prelude::*; 342 /// # let ceed = libceed::Ceed::default_init(); 343 /// let q = 4; 344 /// let b = ceed 345 /// .basis_tensor_H1_Lagrange(2, 1, 3, q, QuadMode::Gauss) 346 /// .unwrap(); 347 /// 348 /// let nqpts = b.num_quadrature_points(); 349 /// assert_eq!(nqpts, q * q, "Incorrect number of quadrature points"); 350 /// ``` 351 pub fn num_quadrature_points(&self) -> usize { 352 let mut Q = 0; 353 unsafe { 354 bind_ceed::CeedBasisGetNumQuadraturePoints(self.ptr, &mut Q); 355 } 356 usize::try_from(Q).unwrap() 357 } 358 } 359 360 // ----------------------------------------------------------------------------- 361