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