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