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 // ElemRestriction 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() -> libceed::Result<()> { 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() -> libceed::Result<()> { 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 // ElemRestriction context wrapper 106 // ----------------------------------------------------------------------------- 107 #[derive(Debug)] 108 pub struct ElemRestriction<'a> { 109 pub(crate) ptr: bind_ceed::CeedElemRestriction, 110 _lifeline: PhantomData<&'a ()>, 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() -> libceed::Result<()> { 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: &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 isize::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 { 200 ptr, 201 _lifeline: PhantomData, 202 }) 203 } 204 205 pub fn create_strided( 206 ceed: &crate::Ceed, 207 nelem: usize, 208 elemsize: usize, 209 ncomp: usize, 210 lsize: usize, 211 strides: [i32; 3], 212 ) -> crate::Result<Self> { 213 let mut ptr = std::ptr::null_mut(); 214 let (nelem, elemsize, ncomp, lsize) = ( 215 i32::try_from(nelem).unwrap(), 216 i32::try_from(elemsize).unwrap(), 217 i32::try_from(ncomp).unwrap(), 218 isize::try_from(lsize).unwrap(), 219 ); 220 let ierr = unsafe { 221 bind_ceed::CeedElemRestrictionCreateStrided( 222 ceed.ptr, 223 nelem, 224 elemsize, 225 ncomp, 226 lsize, 227 strides.as_ptr(), 228 &mut ptr, 229 ) 230 }; 231 ceed.check_error(ierr)?; 232 Ok(Self { 233 ptr, 234 _lifeline: PhantomData, 235 }) 236 } 237 238 // Error handling 239 #[doc(hidden)] 240 fn check_error(&self, ierr: i32) -> crate::Result<i32> { 241 let mut ptr = std::ptr::null_mut(); 242 unsafe { 243 bind_ceed::CeedElemRestrictionGetCeed(self.ptr, &mut ptr); 244 } 245 crate::check_error(ptr, ierr) 246 } 247 248 /// Create an Lvector for an ElemRestriction 249 /// 250 /// ``` 251 /// # use libceed::prelude::*; 252 /// # fn main() -> libceed::Result<()> { 253 /// # let ceed = libceed::Ceed::default_init(); 254 /// let nelem = 3; 255 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 256 /// for i in 0..nelem { 257 /// ind[2 * i + 0] = i as i32; 258 /// ind[2 * i + 1] = (i + 1) as i32; 259 /// } 260 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 261 /// 262 /// let lvector = r.create_lvector()?; 263 /// 264 /// assert_eq!(lvector.length(), nelem + 1, "Incorrect Lvector size"); 265 /// # Ok(()) 266 /// # } 267 /// ``` 268 pub fn create_lvector<'b>(&self) -> crate::Result<Vector<'b>> { 269 let mut ptr_lvector = std::ptr::null_mut(); 270 let null = std::ptr::null_mut() as *mut _; 271 let ierr = 272 unsafe { bind_ceed::CeedElemRestrictionCreateVector(self.ptr, &mut ptr_lvector, null) }; 273 self.check_error(ierr)?; 274 Vector::from_raw(ptr_lvector) 275 } 276 277 /// Create an Evector for an ElemRestriction 278 /// 279 /// ``` 280 /// # use libceed::prelude::*; 281 /// # fn main() -> libceed::Result<()> { 282 /// # let ceed = libceed::Ceed::default_init(); 283 /// let nelem = 3; 284 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 285 /// for i in 0..nelem { 286 /// ind[2 * i + 0] = i as i32; 287 /// ind[2 * i + 1] = (i + 1) as i32; 288 /// } 289 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 290 /// 291 /// let evector = r.create_evector()?; 292 /// 293 /// assert_eq!(evector.length(), nelem * 2, "Incorrect Evector size"); 294 /// # Ok(()) 295 /// # } 296 /// ``` 297 pub fn create_evector<'b>(&self) -> crate::Result<Vector<'b>> { 298 let mut ptr_evector = std::ptr::null_mut(); 299 let null = std::ptr::null_mut() as *mut _; 300 let ierr = 301 unsafe { bind_ceed::CeedElemRestrictionCreateVector(self.ptr, null, &mut ptr_evector) }; 302 self.check_error(ierr)?; 303 Vector::from_raw(ptr_evector) 304 } 305 306 /// Create Vectors for an ElemRestriction 307 /// 308 /// ``` 309 /// # use libceed::prelude::*; 310 /// # fn main() -> libceed::Result<()> { 311 /// # let ceed = libceed::Ceed::default_init(); 312 /// let nelem = 3; 313 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 314 /// for i in 0..nelem { 315 /// ind[2 * i + 0] = i as i32; 316 /// ind[2 * i + 1] = (i + 1) as i32; 317 /// } 318 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 319 /// 320 /// let (lvector, evector) = r.create_vectors()?; 321 /// 322 /// assert_eq!(lvector.length(), nelem + 1, "Incorrect Lvector size"); 323 /// assert_eq!(evector.length(), nelem * 2, "Incorrect Evector size"); 324 /// # Ok(()) 325 /// # } 326 /// ``` 327 pub fn create_vectors<'b, 'c>(&self) -> crate::Result<(Vector<'b>, Vector<'c>)> { 328 let mut ptr_lvector = std::ptr::null_mut(); 329 let mut ptr_evector = std::ptr::null_mut(); 330 let ierr = unsafe { 331 bind_ceed::CeedElemRestrictionCreateVector(self.ptr, &mut ptr_lvector, &mut ptr_evector) 332 }; 333 self.check_error(ierr)?; 334 let lvector = Vector::from_raw(ptr_lvector)?; 335 let evector = Vector::from_raw(ptr_evector)?; 336 Ok((lvector, evector)) 337 } 338 339 /// Restrict an Lvector to an Evector or apply its transpose 340 /// 341 /// # arguments 342 /// 343 /// * `tmode` - Apply restriction or transpose 344 /// * `u` - Input vector (of size `lsize` when `TransposeMode::NoTranspose`) 345 /// * `ru` - Output vector (of shape `[nelem * elemsize]` when 346 /// `TransposeMode::NoTranspose`). Ordering of the Evector is 347 /// decided by the backend. 348 /// 349 /// ``` 350 /// # use libceed::prelude::*; 351 /// # fn main() -> libceed::Result<()> { 352 /// # let ceed = libceed::Ceed::default_init(); 353 /// let nelem = 3; 354 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 355 /// for i in 0..nelem { 356 /// ind[2 * i + 0] = i as i32; 357 /// ind[2 * i + 1] = (i + 1) as i32; 358 /// } 359 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 360 /// 361 /// let x = ceed.vector_from_slice(&[0., 1., 2., 3.])?; 362 /// let mut y = ceed.vector(nelem * 2)?; 363 /// y.set_value(0.0); 364 /// 365 /// r.apply(TransposeMode::NoTranspose, &x, &mut y)?; 366 /// 367 /// for (i, y) in y.view()?.iter().enumerate() { 368 /// assert_eq!( 369 /// *y, 370 /// ((i + 1) / 2) as Scalar, 371 /// "Incorrect value in restricted vector" 372 /// ); 373 /// } 374 /// # Ok(()) 375 /// # } 376 /// ``` 377 pub fn apply(&self, tmode: TransposeMode, u: &Vector, ru: &mut Vector) -> crate::Result<i32> { 378 let tmode = tmode as bind_ceed::CeedTransposeMode; 379 let ierr = unsafe { 380 bind_ceed::CeedElemRestrictionApply( 381 self.ptr, 382 tmode, 383 u.ptr, 384 ru.ptr, 385 bind_ceed::CEED_REQUEST_IMMEDIATE, 386 ) 387 }; 388 self.check_error(ierr) 389 } 390 391 /// Returns the Lvector component stride 392 /// 393 /// ``` 394 /// # use libceed::prelude::*; 395 /// # fn main() -> libceed::Result<()> { 396 /// # let ceed = libceed::Ceed::default_init(); 397 /// let nelem = 3; 398 /// let compstride = 1; 399 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 400 /// for i in 0..nelem { 401 /// ind[2 * i + 0] = i as i32; 402 /// ind[2 * i + 1] = (i + 1) as i32; 403 /// } 404 /// let r = ceed.elem_restriction(nelem, 2, 1, compstride, nelem + 1, MemType::Host, &ind)?; 405 /// 406 /// let c = r.comp_stride(); 407 /// assert_eq!(c, compstride, "Incorrect component stride"); 408 /// # Ok(()) 409 /// # } 410 /// ``` 411 pub fn comp_stride(&self) -> usize { 412 let mut compstride = 0; 413 unsafe { bind_ceed::CeedElemRestrictionGetCompStride(self.ptr, &mut compstride) }; 414 usize::try_from(compstride).unwrap() 415 } 416 417 /// Returns the total number of elements in the range of a ElemRestriction 418 /// 419 /// ``` 420 /// # use libceed::prelude::*; 421 /// # fn main() -> libceed::Result<()> { 422 /// # let ceed = libceed::Ceed::default_init(); 423 /// let nelem = 3; 424 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 425 /// for i in 0..nelem { 426 /// ind[2 * i + 0] = i as i32; 427 /// ind[2 * i + 1] = (i + 1) as i32; 428 /// } 429 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 430 /// 431 /// let n = r.num_elements(); 432 /// assert_eq!(n, nelem, "Incorrect number of elements"); 433 /// # Ok(()) 434 /// # } 435 /// ``` 436 pub fn num_elements(&self) -> usize { 437 let mut numelem = 0; 438 unsafe { bind_ceed::CeedElemRestrictionGetNumElements(self.ptr, &mut numelem) }; 439 usize::try_from(numelem).unwrap() 440 } 441 442 /// Returns the size of elements in the ElemRestriction 443 /// 444 /// ``` 445 /// # use libceed::prelude::*; 446 /// # fn main() -> libceed::Result<()> { 447 /// # let ceed = libceed::Ceed::default_init(); 448 /// let nelem = 3; 449 /// let elem_size = 2; 450 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 451 /// for i in 0..nelem { 452 /// ind[2 * i + 0] = i as i32; 453 /// ind[2 * i + 1] = (i + 1) as i32; 454 /// } 455 /// let r = ceed.elem_restriction(nelem, elem_size, 1, 1, nelem + 1, MemType::Host, &ind)?; 456 /// 457 /// let e = r.elem_size(); 458 /// assert_eq!(e, elem_size, "Incorrect element size"); 459 /// # Ok(()) 460 /// # } 461 /// ``` 462 pub fn elem_size(&self) -> usize { 463 let mut elemsize = 0; 464 unsafe { bind_ceed::CeedElemRestrictionGetElementSize(self.ptr, &mut elemsize) }; 465 usize::try_from(elemsize).unwrap() 466 } 467 468 /// Returns the size of the Lvector for an ElemRestriction 469 /// 470 /// ``` 471 /// # use libceed::prelude::*; 472 /// # fn main() -> libceed::Result<()> { 473 /// # let ceed = libceed::Ceed::default_init(); 474 /// let nelem = 3; 475 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 476 /// for i in 0..nelem { 477 /// ind[2 * i + 0] = i as i32; 478 /// ind[2 * i + 1] = (i + 1) as i32; 479 /// } 480 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 481 /// 482 /// let lsize = r.lvector_size(); 483 /// assert_eq!(lsize, nelem + 1); 484 /// # Ok(()) 485 /// # } 486 /// ``` 487 pub fn lvector_size(&self) -> usize { 488 let mut lsize = 0; 489 unsafe { bind_ceed::CeedElemRestrictionGetLVectorSize(self.ptr, &mut lsize) }; 490 usize::try_from(lsize).unwrap() 491 } 492 493 /// Returns the number of components in the elements of an ElemRestriction 494 /// 495 /// ``` 496 /// # use libceed::prelude::*; 497 /// # fn main() -> libceed::Result<()> { 498 /// # let ceed = libceed::Ceed::default_init(); 499 /// let nelem = 3; 500 /// let ncomp = 42; 501 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 502 /// for i in 0..nelem { 503 /// ind[2 * i + 0] = i as i32; 504 /// ind[2 * i + 1] = (i + 1) as i32; 505 /// } 506 /// let r = ceed.elem_restriction(nelem, 2, 42, 1, ncomp * (nelem + 1), MemType::Host, &ind)?; 507 /// 508 /// let n = r.num_components(); 509 /// assert_eq!(n, ncomp, "Incorrect number of components"); 510 /// # Ok(()) 511 /// # } 512 /// ``` 513 pub fn num_components(&self) -> usize { 514 let mut ncomp = 0; 515 unsafe { bind_ceed::CeedElemRestrictionGetNumComponents(self.ptr, &mut ncomp) }; 516 usize::try_from(ncomp).unwrap() 517 } 518 519 /// Returns the multiplicity of nodes in an ElemRestriction 520 /// 521 /// ``` 522 /// # use libceed::prelude::*; 523 /// # fn main() -> libceed::Result<()> { 524 /// # let ceed = libceed::Ceed::default_init(); 525 /// let nelem = 3; 526 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 527 /// for i in 0..nelem { 528 /// ind[2 * i + 0] = i as i32; 529 /// ind[2 * i + 1] = (i + 1) as i32; 530 /// } 531 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 532 /// 533 /// let mut mult = ceed.vector(nelem + 1)?; 534 /// mult.set_value(0.0); 535 /// 536 /// r.multiplicity(&mut mult)?; 537 /// 538 /// for (i, m) in mult.view()?.iter().enumerate() { 539 /// assert_eq!( 540 /// *m, 541 /// if (i == 0 || i == nelem) { 1. } else { 2. }, 542 /// "Incorrect multiplicity value" 543 /// ); 544 /// } 545 /// # Ok(()) 546 /// # } 547 /// ``` 548 pub fn multiplicity(&self, mult: &mut Vector) -> crate::Result<i32> { 549 let ierr = unsafe { bind_ceed::CeedElemRestrictionGetMultiplicity(self.ptr, mult.ptr) }; 550 self.check_error(ierr) 551 } 552 } 553 554 // ----------------------------------------------------------------------------- 555