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