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 ElemRestriction decomposes elements and groups the degrees of freedom 18 //! (dofs) according to the different elements they belong to. 19 20 use crate::prelude::*; 21 22 // ----------------------------------------------------------------------------- 23 // CeedElemRestriction option 24 // ----------------------------------------------------------------------------- 25 #[derive(Debug)] 26 pub enum ElemRestrictionOpt<'a> { 27 Some(&'a ElemRestriction<'a>), 28 None, 29 } 30 /// Construct a ElemRestrictionOpt reference from a ElemRestriction reference 31 impl<'a> From<&'a ElemRestriction<'_>> for ElemRestrictionOpt<'a> { 32 fn from(restr: &'a ElemRestriction) -> Self { 33 debug_assert!(restr.ptr != unsafe { bind_ceed::CEED_ELEMRESTRICTION_NONE }); 34 Self::Some(restr) 35 } 36 } 37 impl<'a> ElemRestrictionOpt<'a> { 38 /// Transform a Rust libCEED ElemRestrictionOpt into C libCEED 39 /// CeedElemRestriction 40 pub(crate) fn to_raw(self) -> bind_ceed::CeedElemRestriction { 41 match self { 42 Self::Some(restr) => restr.ptr, 43 Self::None => unsafe { bind_ceed::CEED_ELEMRESTRICTION_NONE }, 44 } 45 } 46 } 47 48 // ----------------------------------------------------------------------------- 49 // CeedElemRestriction context wrapper 50 // ----------------------------------------------------------------------------- 51 #[derive(Debug)] 52 pub struct ElemRestriction<'a> { 53 ceed: &'a crate::Ceed, 54 pub(crate) ptr: bind_ceed::CeedElemRestriction, 55 } 56 57 // ----------------------------------------------------------------------------- 58 // Destructor 59 // ----------------------------------------------------------------------------- 60 impl<'a> Drop for ElemRestriction<'a> { 61 fn drop(&mut self) { 62 unsafe { 63 if self.ptr != bind_ceed::CEED_ELEMRESTRICTION_NONE { 64 bind_ceed::CeedElemRestrictionDestroy(&mut self.ptr); 65 } 66 } 67 } 68 } 69 70 // ----------------------------------------------------------------------------- 71 // Display 72 // ----------------------------------------------------------------------------- 73 impl<'a> fmt::Display for ElemRestriction<'a> { 74 /// View an ElemRestriction 75 /// 76 /// ``` 77 /// # use libceed::prelude::*; 78 /// # fn main() -> Result<(), libceed::CeedError> { 79 /// # let ceed = libceed::Ceed::default_init(); 80 /// let nelem = 3; 81 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 82 /// for i in 0..nelem { 83 /// ind[2 * i + 0] = i as i32; 84 /// ind[2 * i + 1] = (i + 1) as i32; 85 /// } 86 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 87 /// println!("{}", r); 88 /// # Ok(()) 89 /// # } 90 /// ``` 91 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 92 let mut ptr = std::ptr::null_mut(); 93 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 94 let cstring = unsafe { 95 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 96 bind_ceed::CeedElemRestrictionView(self.ptr, file); 97 bind_ceed::fclose(file); 98 CString::from_raw(ptr) 99 }; 100 cstring.to_string_lossy().fmt(f) 101 } 102 } 103 104 // ----------------------------------------------------------------------------- 105 // Implementations 106 // ----------------------------------------------------------------------------- 107 impl<'a> ElemRestriction<'a> { 108 // Constructors 109 pub fn create( 110 ceed: &'a crate::Ceed, 111 nelem: usize, 112 elemsize: usize, 113 ncomp: usize, 114 compstride: usize, 115 lsize: usize, 116 mtype: crate::MemType, 117 offsets: &[i32], 118 ) -> crate::Result<Self> { 119 let mut ptr = std::ptr::null_mut(); 120 let (nelem, elemsize, ncomp, compstride, lsize, mtype) = ( 121 i32::try_from(nelem).unwrap(), 122 i32::try_from(elemsize).unwrap(), 123 i32::try_from(ncomp).unwrap(), 124 i32::try_from(compstride).unwrap(), 125 i32::try_from(lsize).unwrap(), 126 mtype as bind_ceed::CeedMemType, 127 ); 128 let ierr = unsafe { 129 bind_ceed::CeedElemRestrictionCreate( 130 ceed.ptr, 131 nelem, 132 elemsize, 133 ncomp, 134 compstride, 135 lsize, 136 mtype, 137 crate::CopyMode::CopyValues as bind_ceed::CeedCopyMode, 138 offsets.as_ptr(), 139 &mut ptr, 140 ) 141 }; 142 ceed.check_error(ierr)?; 143 Ok(Self { ceed, ptr }) 144 } 145 146 pub fn create_strided( 147 ceed: &'a crate::Ceed, 148 nelem: usize, 149 elemsize: usize, 150 ncomp: usize, 151 lsize: usize, 152 strides: [i32; 3], 153 ) -> crate::Result<Self> { 154 let mut ptr = std::ptr::null_mut(); 155 let (nelem, elemsize, ncomp, lsize) = ( 156 i32::try_from(nelem).unwrap(), 157 i32::try_from(elemsize).unwrap(), 158 i32::try_from(ncomp).unwrap(), 159 i32::try_from(lsize).unwrap(), 160 ); 161 let ierr = unsafe { 162 bind_ceed::CeedElemRestrictionCreateStrided( 163 ceed.ptr, 164 nelem, 165 elemsize, 166 ncomp, 167 lsize, 168 strides.as_ptr(), 169 &mut ptr, 170 ) 171 }; 172 ceed.check_error(ierr)?; 173 Ok(Self { ceed, ptr }) 174 } 175 176 /// Create an Lvector for an ElemRestriction 177 /// 178 /// ``` 179 /// # use libceed::prelude::*; 180 /// # fn main() -> Result<(), libceed::CeedError> { 181 /// # let ceed = libceed::Ceed::default_init(); 182 /// let nelem = 3; 183 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 184 /// for i in 0..nelem { 185 /// ind[2 * i + 0] = i as i32; 186 /// ind[2 * i + 1] = (i + 1) as i32; 187 /// } 188 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 189 /// 190 /// let lvector = r.create_lvector()?; 191 /// 192 /// assert_eq!(lvector.length(), nelem + 1, "Incorrect Lvector size"); 193 /// # Ok(()) 194 /// # } 195 /// ``` 196 pub fn create_lvector(&self) -> crate::Result<Vector> { 197 let mut ptr_lvector = std::ptr::null_mut(); 198 let null = std::ptr::null_mut() as *mut _; 199 let ierr = 200 unsafe { bind_ceed::CeedElemRestrictionCreateVector(self.ptr, &mut ptr_lvector, null) }; 201 self.ceed.check_error(ierr)?; 202 Vector::from_raw(self.ceed, ptr_lvector) 203 } 204 205 /// Create an Evector for an ElemRestriction 206 /// 207 /// ``` 208 /// # use libceed::prelude::*; 209 /// # fn main() -> Result<(), libceed::CeedError> { 210 /// # let ceed = libceed::Ceed::default_init(); 211 /// let nelem = 3; 212 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 213 /// for i in 0..nelem { 214 /// ind[2 * i + 0] = i as i32; 215 /// ind[2 * i + 1] = (i + 1) as i32; 216 /// } 217 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 218 /// 219 /// let evector = r.create_evector()?; 220 /// 221 /// assert_eq!(evector.length(), nelem * 2, "Incorrect Evector size"); 222 /// # Ok(()) 223 /// # } 224 /// ``` 225 pub fn create_evector(&self) -> crate::Result<Vector> { 226 let mut ptr_evector = std::ptr::null_mut(); 227 let null = std::ptr::null_mut() as *mut _; 228 let ierr = 229 unsafe { bind_ceed::CeedElemRestrictionCreateVector(self.ptr, null, &mut ptr_evector) }; 230 self.ceed.check_error(ierr)?; 231 Vector::from_raw(self.ceed, ptr_evector) 232 } 233 234 /// Create Vectors for an ElemRestriction 235 /// 236 /// ``` 237 /// # use libceed::prelude::*; 238 /// # fn main() -> Result<(), libceed::CeedError> { 239 /// # let ceed = libceed::Ceed::default_init(); 240 /// let nelem = 3; 241 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 242 /// for i in 0..nelem { 243 /// ind[2 * i + 0] = i as i32; 244 /// ind[2 * i + 1] = (i + 1) as i32; 245 /// } 246 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 247 /// 248 /// let (lvector, evector) = r.create_vectors()?; 249 /// 250 /// assert_eq!(lvector.length(), nelem + 1, "Incorrect Lvector size"); 251 /// assert_eq!(evector.length(), nelem * 2, "Incorrect Evector size"); 252 /// # Ok(()) 253 /// # } 254 /// ``` 255 pub fn create_vectors(&self) -> crate::Result<(Vector, Vector)> { 256 let mut ptr_lvector = std::ptr::null_mut(); 257 let mut ptr_evector = std::ptr::null_mut(); 258 let ierr = unsafe { 259 bind_ceed::CeedElemRestrictionCreateVector(self.ptr, &mut ptr_lvector, &mut ptr_evector) 260 }; 261 self.ceed.check_error(ierr)?; 262 let lvector = Vector::from_raw(self.ceed, ptr_lvector)?; 263 let evector = Vector::from_raw(self.ceed, ptr_evector)?; 264 Ok((lvector, evector)) 265 } 266 267 /// Restrict an Lvector to an Evector or apply its transpose 268 /// 269 /// # arguments 270 /// 271 /// * `tmode` - Apply restriction or transpose 272 /// * `u` - Input vector (of size `lsize` when `TransposeMode::NoTranspose`) 273 /// * `ru` - Output vector (of shape `[nelem * elemsize]` when 274 /// `TransposeMode::NoTranspose`). Ordering of the Evector is 275 /// decided by the backend. 276 /// 277 /// ``` 278 /// # use libceed::prelude::*; 279 /// # fn main() -> Result<(), libceed::CeedError> { 280 /// # let ceed = libceed::Ceed::default_init(); 281 /// let nelem = 3; 282 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 283 /// for i in 0..nelem { 284 /// ind[2 * i + 0] = i as i32; 285 /// ind[2 * i + 1] = (i + 1) as i32; 286 /// } 287 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 288 /// 289 /// let x = ceed.vector_from_slice(&[0., 1., 2., 3.])?; 290 /// let mut y = ceed.vector(nelem * 2)?; 291 /// y.set_value(0.0); 292 /// 293 /// r.apply(TransposeMode::NoTranspose, &x, &mut y)?; 294 /// 295 /// y.view().iter().enumerate().for_each(|(i, arr)| { 296 /// assert_eq!( 297 /// *arr, 298 /// ((i + 1) / 2) as Scalar, 299 /// "Incorrect value in restricted vector" 300 /// ); 301 /// }); 302 /// # Ok(()) 303 /// # } 304 /// ``` 305 pub fn apply(&self, tmode: TransposeMode, u: &Vector, ru: &mut Vector) -> crate::Result<i32> { 306 let tmode = tmode as bind_ceed::CeedTransposeMode; 307 let ierr = unsafe { 308 bind_ceed::CeedElemRestrictionApply( 309 self.ptr, 310 tmode, 311 u.ptr, 312 ru.ptr, 313 bind_ceed::CEED_REQUEST_IMMEDIATE, 314 ) 315 }; 316 self.ceed.check_error(ierr) 317 } 318 319 /// Returns the Lvector component stride 320 /// 321 /// ``` 322 /// # use libceed::prelude::*; 323 /// # fn main() -> Result<(), libceed::CeedError> { 324 /// # let ceed = libceed::Ceed::default_init(); 325 /// let nelem = 3; 326 /// let compstride = 1; 327 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 328 /// for i in 0..nelem { 329 /// ind[2 * i + 0] = i as i32; 330 /// ind[2 * i + 1] = (i + 1) as i32; 331 /// } 332 /// let r = ceed.elem_restriction(nelem, 2, 1, compstride, nelem + 1, MemType::Host, &ind)?; 333 /// 334 /// let c = r.comp_stride(); 335 /// assert_eq!(c, compstride, "Incorrect component stride"); 336 /// # Ok(()) 337 /// # } 338 /// ``` 339 pub fn comp_stride(&self) -> usize { 340 let mut compstride = 0; 341 unsafe { bind_ceed::CeedElemRestrictionGetCompStride(self.ptr, &mut compstride) }; 342 usize::try_from(compstride).unwrap() 343 } 344 345 /// Returns the total number of elements in the range of a ElemRestriction 346 /// 347 /// ``` 348 /// # use libceed::prelude::*; 349 /// # fn main() -> Result<(), libceed::CeedError> { 350 /// # let ceed = libceed::Ceed::default_init(); 351 /// let nelem = 3; 352 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 353 /// for i in 0..nelem { 354 /// ind[2 * i + 0] = i as i32; 355 /// ind[2 * i + 1] = (i + 1) as i32; 356 /// } 357 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 358 /// 359 /// let n = r.num_elements(); 360 /// assert_eq!(n, nelem, "Incorrect number of elements"); 361 /// # Ok(()) 362 /// # } 363 /// ``` 364 pub fn num_elements(&self) -> usize { 365 let mut numelem = 0; 366 unsafe { bind_ceed::CeedElemRestrictionGetNumElements(self.ptr, &mut numelem) }; 367 usize::try_from(numelem).unwrap() 368 } 369 370 /// Returns the size of elements in the ElemRestriction 371 /// 372 /// ``` 373 /// # use libceed::prelude::*; 374 /// # fn main() -> Result<(), libceed::CeedError> { 375 /// # let ceed = libceed::Ceed::default_init(); 376 /// let nelem = 3; 377 /// let elem_size = 2; 378 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 379 /// for i in 0..nelem { 380 /// ind[2 * i + 0] = i as i32; 381 /// ind[2 * i + 1] = (i + 1) as i32; 382 /// } 383 /// let r = ceed.elem_restriction(nelem, elem_size, 1, 1, nelem + 1, MemType::Host, &ind)?; 384 /// 385 /// let e = r.elem_size(); 386 /// assert_eq!(e, elem_size, "Incorrect element size"); 387 /// # Ok(()) 388 /// # } 389 /// ``` 390 pub fn elem_size(&self) -> usize { 391 let mut elemsize = 0; 392 unsafe { bind_ceed::CeedElemRestrictionGetElementSize(self.ptr, &mut elemsize) }; 393 usize::try_from(elemsize).unwrap() 394 } 395 396 /// Returns the size of the Lvector for an ElemRestriction 397 /// 398 /// ``` 399 /// # use libceed::prelude::*; 400 /// # fn main() -> Result<(), libceed::CeedError> { 401 /// # let ceed = libceed::Ceed::default_init(); 402 /// let nelem = 3; 403 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 404 /// for i in 0..nelem { 405 /// ind[2 * i + 0] = i as i32; 406 /// ind[2 * i + 1] = (i + 1) as i32; 407 /// } 408 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 409 /// 410 /// let lsize = r.lvector_size(); 411 /// assert_eq!(lsize, nelem + 1); 412 /// # Ok(()) 413 /// # } 414 /// ``` 415 pub fn lvector_size(&self) -> usize { 416 let mut lsize = 0; 417 unsafe { bind_ceed::CeedElemRestrictionGetLVectorSize(self.ptr, &mut lsize) }; 418 usize::try_from(lsize).unwrap() 419 } 420 421 /// Returns the number of components in the elements of an ElemRestriction 422 /// 423 /// ``` 424 /// # use libceed::prelude::*; 425 /// # fn main() -> Result<(), libceed::CeedError> { 426 /// # let ceed = libceed::Ceed::default_init(); 427 /// let nelem = 3; 428 /// let ncomp = 42; 429 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 430 /// for i in 0..nelem { 431 /// ind[2 * i + 0] = i as i32; 432 /// ind[2 * i + 1] = (i + 1) as i32; 433 /// } 434 /// let r = ceed.elem_restriction(nelem, 2, 42, 1, ncomp * (nelem + 1), MemType::Host, &ind)?; 435 /// 436 /// let n = r.num_components(); 437 /// assert_eq!(n, ncomp, "Incorrect number of components"); 438 /// # Ok(()) 439 /// # } 440 /// ``` 441 pub fn num_components(&self) -> usize { 442 let mut ncomp = 0; 443 unsafe { bind_ceed::CeedElemRestrictionGetNumComponents(self.ptr, &mut ncomp) }; 444 usize::try_from(ncomp).unwrap() 445 } 446 447 /// Returns the multiplicity of nodes in an ElemRestriction 448 /// 449 /// ``` 450 /// # use libceed::prelude::*; 451 /// # fn main() -> Result<(), libceed::CeedError> { 452 /// # let ceed = libceed::Ceed::default_init(); 453 /// let nelem = 3; 454 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 455 /// for i in 0..nelem { 456 /// ind[2 * i + 0] = i as i32; 457 /// ind[2 * i + 1] = (i + 1) as i32; 458 /// } 459 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 460 /// 461 /// let mut mult = ceed.vector(nelem + 1)?; 462 /// mult.set_value(0.0); 463 /// 464 /// r.multiplicity(&mut mult)?; 465 /// 466 /// mult.view().iter().enumerate().for_each(|(i, arr)| { 467 /// assert_eq!( 468 /// if (i == 0 || i == nelem) { 1. } else { 2. }, 469 /// *arr, 470 /// "Incorrect multiplicity array" 471 /// ); 472 /// }); 473 /// # Ok(()) 474 /// # } 475 /// ``` 476 pub fn multiplicity(&self, mult: &mut Vector) -> crate::Result<i32> { 477 let ierr = unsafe { bind_ceed::CeedElemRestrictionGetMultiplicity(self.ptr, mult.ptr) }; 478 self.ceed.check_error(ierr) 479 } 480 } 481 482 // ----------------------------------------------------------------------------- 483