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: &[f64], 107 grad1d: &[f64], 108 qref1d: &[f64], 109 qweight1d: &[f64], 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: &[f64], 166 grad: &[f64], 167 qref: &[f64], 168 qweight: &[f64], 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_eq!(*v, true_value, "Incorrect basis application"); 257 /// }); 258 /// ``` 259 pub fn apply( 260 &self, 261 nelem: usize, 262 tmode: TransposeMode, 263 emode: EvalMode, 264 u: &Vector, 265 v: &mut Vector, 266 ) -> crate::Result<i32> { 267 let (nelem, tmode, emode) = ( 268 i32::try_from(nelem).unwrap(), 269 tmode as bind_ceed::CeedTransposeMode, 270 emode as bind_ceed::CeedEvalMode, 271 ); 272 let ierr = 273 unsafe { bind_ceed::CeedBasisApply(self.ptr, nelem, tmode, emode, u.ptr, v.ptr) }; 274 self.ceed.check_error(ierr) 275 } 276 277 /// Returns the dimension for given CeedBasis 278 /// 279 /// ``` 280 /// # use libceed::prelude::*; 281 /// # let ceed = libceed::Ceed::default_init(); 282 /// let dim = 2; 283 /// let b = ceed 284 /// .basis_tensor_H1_Lagrange(dim, 1, 3, 4, QuadMode::Gauss) 285 /// .unwrap(); 286 /// 287 /// let d = b.dimension(); 288 /// assert_eq!(d, dim, "Incorrect dimension"); 289 /// ``` 290 pub fn dimension(&self) -> usize { 291 let mut dim = 0; 292 unsafe { bind_ceed::CeedBasisGetDimension(self.ptr, &mut dim) }; 293 usize::try_from(dim).unwrap() 294 } 295 296 /// Returns number of components for given CeedBasis 297 /// 298 /// ``` 299 /// # use libceed::prelude::*; 300 /// # let ceed = libceed::Ceed::default_init(); 301 /// let ncomp = 2; 302 /// let b = ceed 303 /// .basis_tensor_H1_Lagrange(1, ncomp, 3, 4, QuadMode::Gauss) 304 /// .unwrap(); 305 /// 306 /// let n = b.num_components(); 307 /// assert_eq!(n, ncomp, "Incorrect number of components"); 308 /// ``` 309 pub fn num_components(&self) -> usize { 310 let mut ncomp = 0; 311 unsafe { bind_ceed::CeedBasisGetNumComponents(self.ptr, &mut ncomp) }; 312 usize::try_from(ncomp).unwrap() 313 } 314 315 /// Returns total number of nodes (in dim dimensions) of a CeedBasis 316 /// 317 /// ``` 318 /// # use libceed::prelude::*; 319 /// # let ceed = libceed::Ceed::default_init(); 320 /// let p = 3; 321 /// let b = ceed 322 /// .basis_tensor_H1_Lagrange(2, 1, p, 4, QuadMode::Gauss) 323 /// .unwrap(); 324 /// 325 /// let nnodes = b.num_nodes(); 326 /// assert_eq!(nnodes, p * p, "Incorrect number of nodes"); 327 /// ``` 328 pub fn num_nodes(&self) -> usize { 329 let mut nnodes = 0; 330 unsafe { bind_ceed::CeedBasisGetNumNodes(self.ptr, &mut nnodes) }; 331 usize::try_from(nnodes).unwrap() 332 } 333 334 /// Returns total number of quadrature points (in dim dimensions) of a 335 /// CeedBasis 336 /// 337 /// ``` 338 /// # use libceed::prelude::*; 339 /// # let ceed = libceed::Ceed::default_init(); 340 /// let q = 4; 341 /// let b = ceed 342 /// .basis_tensor_H1_Lagrange(2, 1, 3, q, QuadMode::Gauss) 343 /// .unwrap(); 344 /// 345 /// let nqpts = b.num_quadrature_points(); 346 /// assert_eq!(nqpts, q * q, "Incorrect number of quadrature points"); 347 /// ``` 348 pub fn num_quadrature_points(&self) -> usize { 349 let mut Q = 0; 350 unsafe { 351 bind_ceed::CeedBasisGetNumQuadraturePoints(self.ptr, &mut Q); 352 } 353 usize::try_from(Q).unwrap() 354 } 355 } 356 357 // ----------------------------------------------------------------------------- 358