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 // Basis 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() -> libceed::Result<()> { 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() -> libceed::Result<()> { 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 // Basis context wrapper 93 // ----------------------------------------------------------------------------- 94 #[derive(Debug)] 95 pub struct Basis<'a> { 96 pub(crate) ptr: bind_ceed::CeedBasis, 97 _lifeline: PhantomData<&'a ()>, 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() -> libceed::Result<()> { 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: &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 { 180 ptr, 181 _lifeline: PhantomData, 182 }) 183 } 184 185 pub fn create_tensor_H1_Lagrange( 186 ceed: &crate::Ceed, 187 dim: usize, 188 ncomp: usize, 189 P: usize, 190 Q: usize, 191 qmode: crate::QuadMode, 192 ) -> crate::Result<Self> { 193 let mut ptr = std::ptr::null_mut(); 194 let (dim, ncomp, P, Q, qmode) = ( 195 i32::try_from(dim).unwrap(), 196 i32::try_from(ncomp).unwrap(), 197 i32::try_from(P).unwrap(), 198 i32::try_from(Q).unwrap(), 199 qmode as bind_ceed::CeedQuadMode, 200 ); 201 let ierr = unsafe { 202 bind_ceed::CeedBasisCreateTensorH1Lagrange(ceed.ptr, dim, ncomp, P, Q, qmode, &mut ptr) 203 }; 204 ceed.check_error(ierr)?; 205 Ok(Self { 206 ptr, 207 _lifeline: PhantomData, 208 }) 209 } 210 211 pub fn create_H1( 212 ceed: &crate::Ceed, 213 topo: crate::ElemTopology, 214 ncomp: usize, 215 nnodes: usize, 216 nqpts: usize, 217 interp: &[crate::Scalar], 218 grad: &[crate::Scalar], 219 qref: &[crate::Scalar], 220 qweight: &[crate::Scalar], 221 ) -> crate::Result<Self> { 222 let mut ptr = std::ptr::null_mut(); 223 let (topo, ncomp, nnodes, nqpts) = ( 224 topo as bind_ceed::CeedElemTopology, 225 i32::try_from(ncomp).unwrap(), 226 i32::try_from(nnodes).unwrap(), 227 i32::try_from(nqpts).unwrap(), 228 ); 229 let ierr = unsafe { 230 bind_ceed::CeedBasisCreateH1( 231 ceed.ptr, 232 topo, 233 ncomp, 234 nnodes, 235 nqpts, 236 interp.as_ptr(), 237 grad.as_ptr(), 238 qref.as_ptr(), 239 qweight.as_ptr(), 240 &mut ptr, 241 ) 242 }; 243 ceed.check_error(ierr)?; 244 Ok(Self { 245 ptr, 246 _lifeline: PhantomData, 247 }) 248 } 249 250 // Error handling 251 #[doc(hidden)] 252 fn check_error(&self, ierr: i32) -> crate::Result<i32> { 253 let mut ptr = std::ptr::null_mut(); 254 unsafe { 255 bind_ceed::CeedBasisGetCeed(self.ptr, &mut ptr); 256 } 257 crate::check_error(ptr, ierr) 258 } 259 260 /// Apply basis evaluation from nodes to quadrature points or vice versa 261 /// 262 /// * `nelem` - The number of elements to apply the basis evaluation to 263 /// * `tmode` - `TrasposeMode::NoTranspose` to evaluate from nodes to 264 /// quadrature points, `TransposeMode::Transpose` to apply the 265 /// transpose, mapping from quadrature points to nodes 266 /// * `emode` - `EvalMode::None` to use values directly, `EvalMode::Interp` 267 /// to use interpolated values, `EvalMode::Grad` to use 268 /// gradients, `EvalMode::Weight` to use quadrature weights 269 /// * `u` - Input Vector 270 /// * `v` - Output Vector 271 /// 272 /// ``` 273 /// # use libceed::prelude::*; 274 /// # fn main() -> libceed::Result<()> { 275 /// # let ceed = libceed::Ceed::default_init(); 276 /// const Q: usize = 6; 277 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, Q, Q, QuadMode::GaussLobatto)?; 278 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, Q, QuadMode::Gauss)?; 279 /// 280 /// let x_corners = ceed.vector_from_slice(&[-1., 1.])?; 281 /// let mut x_qpts = ceed.vector(Q)?; 282 /// let mut x_nodes = ceed.vector(Q)?; 283 /// bx.apply( 284 /// 1, 285 /// TransposeMode::NoTranspose, 286 /// EvalMode::Interp, 287 /// &x_corners, 288 /// &mut x_nodes, 289 /// )?; 290 /// bu.apply( 291 /// 1, 292 /// TransposeMode::NoTranspose, 293 /// EvalMode::Interp, 294 /// &x_nodes, 295 /// &mut x_qpts, 296 /// )?; 297 /// 298 /// // Create function x^3 + 1 on Gauss Lobatto points 299 /// let mut u_arr = [0.; Q]; 300 /// u_arr 301 /// .iter_mut() 302 /// .zip(x_nodes.view()?.iter()) 303 /// .for_each(|(u, x)| *u = x * x * x + 1.); 304 /// let u = ceed.vector_from_slice(&u_arr)?; 305 /// 306 /// // Map function to Gauss points 307 /// let mut v = ceed.vector(Q)?; 308 /// v.set_value(0.); 309 /// bu.apply(1, TransposeMode::NoTranspose, EvalMode::Interp, &u, &mut v)?; 310 /// 311 /// // Verify results 312 /// v.view()? 313 /// .iter() 314 /// .zip(x_qpts.view()?.iter()) 315 /// .for_each(|(v, x)| { 316 /// let true_value = x * x * x + 1.; 317 /// assert!( 318 /// (*v - true_value).abs() < 10.0 * libceed::EPSILON, 319 /// "Incorrect basis application" 320 /// ); 321 /// }); 322 /// # Ok(()) 323 /// # } 324 /// ``` 325 pub fn apply( 326 &self, 327 nelem: usize, 328 tmode: TransposeMode, 329 emode: EvalMode, 330 u: &Vector, 331 v: &mut Vector, 332 ) -> crate::Result<i32> { 333 let (nelem, tmode, emode) = ( 334 i32::try_from(nelem).unwrap(), 335 tmode as bind_ceed::CeedTransposeMode, 336 emode as bind_ceed::CeedEvalMode, 337 ); 338 let ierr = 339 unsafe { bind_ceed::CeedBasisApply(self.ptr, nelem, tmode, emode, u.ptr, v.ptr) }; 340 self.check_error(ierr) 341 } 342 343 /// Returns the dimension for given Basis 344 /// 345 /// ``` 346 /// # use libceed::prelude::*; 347 /// # fn main() -> libceed::Result<()> { 348 /// # let ceed = libceed::Ceed::default_init(); 349 /// let dim = 2; 350 /// let b = ceed.basis_tensor_H1_Lagrange(dim, 1, 3, 4, QuadMode::Gauss)?; 351 /// 352 /// let d = b.dimension(); 353 /// assert_eq!(d, dim, "Incorrect dimension"); 354 /// # Ok(()) 355 /// # } 356 /// ``` 357 pub fn dimension(&self) -> usize { 358 let mut dim = 0; 359 unsafe { bind_ceed::CeedBasisGetDimension(self.ptr, &mut dim) }; 360 usize::try_from(dim).unwrap() 361 } 362 363 /// Returns number of components for given Basis 364 /// 365 /// ``` 366 /// # use libceed::prelude::*; 367 /// # fn main() -> libceed::Result<()> { 368 /// # let ceed = libceed::Ceed::default_init(); 369 /// let ncomp = 2; 370 /// let b = ceed.basis_tensor_H1_Lagrange(1, ncomp, 3, 4, QuadMode::Gauss)?; 371 /// 372 /// let n = b.num_components(); 373 /// assert_eq!(n, ncomp, "Incorrect number of components"); 374 /// # Ok(()) 375 /// # } 376 /// ``` 377 pub fn num_components(&self) -> usize { 378 let mut ncomp = 0; 379 unsafe { bind_ceed::CeedBasisGetNumComponents(self.ptr, &mut ncomp) }; 380 usize::try_from(ncomp).unwrap() 381 } 382 383 /// Returns total number of nodes (in dim dimensions) of a Basis 384 /// 385 /// ``` 386 /// # use libceed::prelude::*; 387 /// # fn main() -> libceed::Result<()> { 388 /// # let ceed = libceed::Ceed::default_init(); 389 /// let p = 3; 390 /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, p, 4, QuadMode::Gauss)?; 391 /// 392 /// let nnodes = b.num_nodes(); 393 /// assert_eq!(nnodes, p * p, "Incorrect number of nodes"); 394 /// # Ok(()) 395 /// # } 396 /// ``` 397 pub fn num_nodes(&self) -> usize { 398 let mut nnodes = 0; 399 unsafe { bind_ceed::CeedBasisGetNumNodes(self.ptr, &mut nnodes) }; 400 usize::try_from(nnodes).unwrap() 401 } 402 403 /// Returns total number of quadrature points (in dim dimensions) of a 404 /// Basis 405 /// 406 /// ``` 407 /// # use libceed::prelude::*; 408 /// # fn main() -> libceed::Result<()> { 409 /// # let ceed = libceed::Ceed::default_init(); 410 /// let q = 4; 411 /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, 3, q, QuadMode::Gauss)?; 412 /// 413 /// let nqpts = b.num_quadrature_points(); 414 /// assert_eq!(nqpts, q * q, "Incorrect number of quadrature points"); 415 /// # Ok(()) 416 /// # } 417 /// ``` 418 pub fn num_quadrature_points(&self) -> usize { 419 let mut Q = 0; 420 unsafe { 421 bind_ceed::CeedBasisGetNumQuadraturePoints(self.ptr, &mut Q); 422 } 423 usize::try_from(Q).unwrap() 424 } 425 } 426 427 // ----------------------------------------------------------------------------- 428