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 // Fenced `rust` code blocks included from README.md are executed as part of doctests. 18 #![doc = include_str!("../README.md")] 19 // ----------------------------------------------------------------------------- 20 // Exceptions 21 // ----------------------------------------------------------------------------- 22 #![allow(non_snake_case)] 23 24 // ----------------------------------------------------------------------------- 25 // Crate prelude 26 // ----------------------------------------------------------------------------- 27 use crate::prelude::*; 28 use std::sync::Once; 29 30 pub mod prelude { 31 pub use crate::{ 32 basis::{self, Basis, BasisOpt}, 33 elem_restriction::{self, ElemRestriction, ElemRestrictionOpt}, 34 operator::{self, CompositeOperator, Operator}, 35 qfunction::{ 36 self, QFunction, QFunctionByName, QFunctionInputs, QFunctionOpt, QFunctionOutputs, 37 }, 38 vector::{self, Vector, VectorOpt}, 39 ElemTopology, EvalMode, MemType, NormType, QuadMode, Scalar, TransposeMode, 40 CEED_STRIDES_BACKEND, EPSILON, MAX_QFUNCTION_FIELDS, 41 }; 42 pub(crate) use libceed_sys::bind_ceed; 43 pub(crate) use std::convert::TryFrom; 44 pub(crate) use std::ffi::CString; 45 pub(crate) use std::fmt; 46 } 47 48 // ----------------------------------------------------------------------------- 49 // Modules 50 // ----------------------------------------------------------------------------- 51 pub mod basis; 52 pub mod elem_restriction; 53 pub mod operator; 54 pub mod qfunction; 55 pub mod vector; 56 57 // ----------------------------------------------------------------------------- 58 // Typedef for scalar 59 // ----------------------------------------------------------------------------- 60 pub type Scalar = bind_ceed::CeedScalar; 61 62 // ----------------------------------------------------------------------------- 63 // Constants for library 64 // ----------------------------------------------------------------------------- 65 const MAX_BUFFER_LENGTH: u64 = 4096; 66 pub const MAX_QFUNCTION_FIELDS: usize = 16; 67 pub const CEED_STRIDES_BACKEND: [i32; 3] = [0; 3]; 68 pub const EPSILON: crate::Scalar = bind_ceed::CEED_EPSILON as crate::Scalar; 69 70 // ----------------------------------------------------------------------------- 71 // Enums for libCEED 72 // ----------------------------------------------------------------------------- 73 #[derive(Clone, Copy, PartialEq, Eq)] 74 /// Many Ceed interfaces take or return pointers to memory. This enum is used to 75 /// specify where the memory being provided or requested must reside. 76 pub enum MemType { 77 Host = bind_ceed::CeedMemType_CEED_MEM_HOST as isize, 78 Device = bind_ceed::CeedMemType_CEED_MEM_DEVICE as isize, 79 } 80 81 #[derive(Clone, Copy, PartialEq, Eq)] 82 // OwnPointer will not be used by user but is included for internal use 83 #[allow(dead_code)] 84 /// Conveys ownership status of arrays passed to Ceed interfaces. 85 pub(crate) enum CopyMode { 86 CopyValues = bind_ceed::CeedCopyMode_CEED_COPY_VALUES as isize, 87 UsePointer = bind_ceed::CeedCopyMode_CEED_USE_POINTER as isize, 88 OwnPointer = bind_ceed::CeedCopyMode_CEED_OWN_POINTER as isize, 89 } 90 91 #[derive(Clone, Copy, PartialEq, Eq)] 92 /// Denotes type of vector norm to be computed 93 pub enum NormType { 94 One = bind_ceed::CeedNormType_CEED_NORM_1 as isize, 95 Two = bind_ceed::CeedNormType_CEED_NORM_2 as isize, 96 Max = bind_ceed::CeedNormType_CEED_NORM_MAX as isize, 97 } 98 99 #[derive(Clone, Copy, PartialEq, Eq)] 100 /// Denotes whether a linear transformation or its transpose should be applied 101 pub enum TransposeMode { 102 NoTranspose = bind_ceed::CeedTransposeMode_CEED_NOTRANSPOSE as isize, 103 Transpose = bind_ceed::CeedTransposeMode_CEED_TRANSPOSE as isize, 104 } 105 106 #[derive(Clone, Copy, PartialEq, Eq)] 107 /// Type of quadrature; also used for location of nodes 108 pub enum QuadMode { 109 Gauss = bind_ceed::CeedQuadMode_CEED_GAUSS as isize, 110 GaussLobatto = bind_ceed::CeedQuadMode_CEED_GAUSS_LOBATTO as isize, 111 } 112 113 #[derive(Clone, Copy, PartialEq, Eq)] 114 /// Type of basis shape to create non-tensor H1 element basis 115 pub enum ElemTopology { 116 Line = bind_ceed::CeedElemTopology_CEED_LINE as isize, 117 Triangle = bind_ceed::CeedElemTopology_CEED_TRIANGLE as isize, 118 Quad = bind_ceed::CeedElemTopology_CEED_QUAD as isize, 119 Tet = bind_ceed::CeedElemTopology_CEED_TET as isize, 120 Pyramid = bind_ceed::CeedElemTopology_CEED_PYRAMID as isize, 121 Prism = bind_ceed::CeedElemTopology_CEED_PRISM as isize, 122 Hex = bind_ceed::CeedElemTopology_CEED_HEX as isize, 123 } 124 125 #[derive(Clone, Copy, PartialEq, Eq)] 126 /// Basis evaluation mode 127 pub enum EvalMode { 128 None = bind_ceed::CeedEvalMode_CEED_EVAL_NONE as isize, 129 Interp = bind_ceed::CeedEvalMode_CEED_EVAL_INTERP as isize, 130 Grad = bind_ceed::CeedEvalMode_CEED_EVAL_GRAD as isize, 131 Div = bind_ceed::CeedEvalMode_CEED_EVAL_DIV as isize, 132 Curl = bind_ceed::CeedEvalMode_CEED_EVAL_CURL as isize, 133 Weight = bind_ceed::CeedEvalMode_CEED_EVAL_WEIGHT as isize, 134 } 135 136 // ----------------------------------------------------------------------------- 137 // Ceed error 138 // ----------------------------------------------------------------------------- 139 type Result<T> = std::result::Result<T, CeedError>; 140 141 #[derive(Debug)] 142 pub struct CeedError { 143 pub message: String, 144 } 145 146 impl fmt::Display for CeedError { 147 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 148 write!(f, "{}", self.message) 149 } 150 } 151 152 // ----------------------------------------------------------------------------- 153 // Ceed error handler 154 // ----------------------------------------------------------------------------- 155 pub enum CeedErrorHandler { 156 ErrorAbort, 157 ErrorExit, 158 ErrorReturn, 159 ErrorStore, 160 } 161 162 // ----------------------------------------------------------------------------- 163 // Ceed context wrapper 164 // ----------------------------------------------------------------------------- 165 /// A Ceed is a library context representing control of a logical hardware 166 /// resource. 167 #[derive(Debug)] 168 pub struct Ceed { 169 ptr: bind_ceed::Ceed, 170 } 171 172 // ----------------------------------------------------------------------------- 173 // Destructor 174 // ----------------------------------------------------------------------------- 175 impl Drop for Ceed { 176 fn drop(&mut self) { 177 unsafe { 178 bind_ceed::CeedDestroy(&mut self.ptr); 179 } 180 } 181 } 182 183 // ----------------------------------------------------------------------------- 184 // Display 185 // ----------------------------------------------------------------------------- 186 impl fmt::Display for Ceed { 187 /// View a Ceed 188 /// 189 /// ``` 190 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 191 /// println!("{}", ceed); 192 /// ``` 193 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 194 let mut ptr = std::ptr::null_mut(); 195 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 196 let cstring = unsafe { 197 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 198 bind_ceed::CeedView(self.ptr, file); 199 bind_ceed::fclose(file); 200 CString::from_raw(ptr) 201 }; 202 cstring.to_string_lossy().fmt(f) 203 } 204 } 205 206 static REGISTER: Once = Once::new(); 207 208 // ----------------------------------------------------------------------------- 209 // Object constructors 210 // ----------------------------------------------------------------------------- 211 impl Ceed { 212 /// Returns a Ceed context initialized with the specified resource 213 /// 214 /// # arguments 215 /// 216 /// * `resource` - Resource to use, e.g., "/cpu/self" 217 /// 218 /// ``` 219 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 220 /// ``` 221 pub fn init(resource: &str) -> Self { 222 Ceed::init_with_error_handler(resource, CeedErrorHandler::ErrorStore) 223 } 224 225 /// Returns a Ceed context initialized with the specified resource 226 /// 227 /// # arguments 228 /// 229 /// * `resource` - Resource to use, e.g., "/cpu/self" 230 /// 231 /// ``` 232 /// let ceed = libceed::Ceed::init_with_error_handler( 233 /// "/cpu/self/ref/serial", 234 /// libceed::CeedErrorHandler::ErrorAbort, 235 /// ); 236 /// ``` 237 pub fn init_with_error_handler(resource: &str, handler: CeedErrorHandler) -> Self { 238 REGISTER.call_once(|| unsafe { 239 bind_ceed::CeedRegisterAll(); 240 bind_ceed::CeedQFunctionRegisterAll(); 241 }); 242 243 // Convert to C string 244 let c_resource = CString::new(resource).expect("CString::new failed"); 245 246 // Get error handler pointer 247 let eh = match handler { 248 CeedErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort, 249 CeedErrorHandler::ErrorExit => bind_ceed::CeedErrorExit, 250 CeedErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn, 251 CeedErrorHandler::ErrorStore => bind_ceed::CeedErrorStore, 252 }; 253 254 // Call to libCEED 255 let mut ptr = std::ptr::null_mut(); 256 let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr() as *const i8, &mut ptr) }; 257 if ierr != 0 { 258 panic!("Error initializing backend resource: {}", resource) 259 } 260 ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) }; 261 let ceed = Ceed { ptr }; 262 ceed.check_error(ierr).unwrap(); 263 ceed 264 } 265 266 /// Default initializer for testing 267 #[doc(hidden)] 268 pub fn default_init() -> Self { 269 // Convert to C string 270 let resource = "/cpu/self/ref/serial"; 271 crate::Ceed::init(resource) 272 } 273 274 /// Internal error checker 275 #[doc(hidden)] 276 fn check_error(&self, ierr: i32) -> Result<i32> { 277 // Return early if code is clean 278 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 279 return Ok(ierr); 280 } 281 // Retrieve error message 282 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 283 let c_str = unsafe { 284 bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr); 285 std::ffi::CStr::from_ptr(ptr) 286 }; 287 let message = c_str.to_string_lossy().to_string(); 288 // Panic if negative code, otherwise return error 289 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 290 panic!("{}", message); 291 } 292 Err(CeedError { message }) 293 } 294 295 /// Returns full resource name for a Ceed context 296 /// 297 /// ``` 298 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 299 /// let resource = ceed.resource(); 300 /// 301 /// assert_eq!(resource, "/cpu/self/ref/serial".to_string()) 302 /// ``` 303 pub fn resource(&self) -> String { 304 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 305 let c_str = unsafe { 306 bind_ceed::CeedGetResource(self.ptr, &mut ptr); 307 std::ffi::CStr::from_ptr(ptr) 308 }; 309 c_str.to_string_lossy().to_string() 310 } 311 312 /// Returns a CeedVector of the specified length (does not allocate memory) 313 /// 314 /// # arguments 315 /// 316 /// * `n` - Length of vector 317 /// 318 /// ``` 319 /// # use libceed::prelude::*; 320 /// # let ceed = libceed::Ceed::default_init(); 321 /// let vec = ceed.vector(10).unwrap(); 322 /// ``` 323 pub fn vector(&self, n: usize) -> Result<Vector> { 324 Vector::create(self, n) 325 } 326 327 /// Create a Vector initialized with the data (copied) from a slice 328 /// 329 /// # arguments 330 /// 331 /// * `slice` - Slice containing data 332 /// 333 /// ``` 334 /// # use libceed::prelude::*; 335 /// # let ceed = libceed::Ceed::default_init(); 336 /// let vec = ceed.vector_from_slice(&[1., 2., 3.]).unwrap(); 337 /// assert_eq!(vec.length(), 3); 338 /// ``` 339 pub fn vector_from_slice(&self, slice: &[crate::Scalar]) -> Result<Vector> { 340 Vector::from_slice(self, slice) 341 } 342 343 /// Returns a ElemRestriction 344 /// 345 /// # arguments 346 /// 347 /// * `nelem` - Number of elements described in the offsets array 348 /// * `elemsize` - Size (number of "nodes") per element 349 /// * `ncomp` - Number of field components per interpolation node (1 350 /// for scalar fields) 351 /// * `compstride` - Stride between components for the same Lvector "node". 352 /// Data for node `i`, component `j`, element `k` can be 353 /// found in the Lvector at index 354 /// `offsets[i + k*elemsize] + j*compstride`. 355 /// * `lsize` - The size of the Lvector. This vector may be larger 356 /// than the elements and fields given by this 357 /// restriction. 358 /// * `mtype` - Memory type of the offsets array, see CeedMemType 359 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 360 /// ordered list of the offsets (into the input CeedVector) 361 /// for the unknowns corresponding to element `i`, where 362 /// `0 <= i < nelem`. All offsets must be in the range 363 /// `[0, lsize - 1]`. 364 /// 365 /// ``` 366 /// # use libceed::prelude::*; 367 /// # let ceed = libceed::Ceed::default_init(); 368 /// let nelem = 3; 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, 2, 1, 1, nelem + 1, MemType::Host, &ind) 376 /// .unwrap(); 377 /// ``` 378 pub fn elem_restriction( 379 &self, 380 nelem: usize, 381 elemsize: usize, 382 ncomp: usize, 383 compstride: usize, 384 lsize: usize, 385 mtype: MemType, 386 offsets: &[i32], 387 ) -> Result<ElemRestriction> { 388 ElemRestriction::create( 389 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, 390 ) 391 } 392 393 /// Returns a ElemRestriction 394 /// 395 /// # arguments 396 /// 397 /// * `nelem` - Number of elements described in the offsets array 398 /// * `elemsize` - Size (number of "nodes") per element 399 /// * `ncomp` - Number of field components per interpolation node (1 400 /// for scalar fields) 401 /// * `compstride` - Stride between components for the same Lvector "node". 402 /// Data for node `i`, component `j`, element `k` can be 403 /// found in the Lvector at index 404 /// `offsets[i + k*elemsize] + j*compstride`. 405 /// * `lsize` - The size of the Lvector. This vector may be larger 406 /// than the elements and fields given by this restriction. 407 /// * `strides` - Array for strides between `[nodes, components, elements]`. 408 /// Data for node `i`, component `j`, element `k` can be 409 /// found in the Lvector at index 410 /// `i*strides[0] + j*strides[1] + k*strides[2]`. 411 /// CEED_STRIDES_BACKEND may be used with vectors created 412 /// by a Ceed backend. 413 /// 414 /// ``` 415 /// # use libceed::prelude::*; 416 /// # let ceed = libceed::Ceed::default_init(); 417 /// let nelem = 3; 418 /// let strides: [i32; 3] = [1, 2, 2]; 419 /// let r = ceed 420 /// .strided_elem_restriction(nelem, 2, 1, nelem * 2, strides) 421 /// .unwrap(); 422 /// ``` 423 pub fn strided_elem_restriction( 424 &self, 425 nelem: usize, 426 elemsize: usize, 427 ncomp: usize, 428 lsize: usize, 429 strides: [i32; 3], 430 ) -> Result<ElemRestriction> { 431 ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides) 432 } 433 434 /// Returns a tensor-product basis 435 /// 436 /// # arguments 437 /// 438 /// * `dim` - Topological dimension of element 439 /// * `ncomp` - Number of field components (1 for scalar fields) 440 /// * `P1d` - Number of Gauss-Lobatto nodes in one dimension. The 441 /// polynomial degree of the resulting `Q_k` element is 442 /// `k=P-1`. 443 /// * `Q1d` - Number of quadrature points in one dimension 444 /// * `interp1d` - Row-major `(Q1d * P1d)` matrix expressing the values of 445 /// nodal basis functions at quadrature points 446 /// * `grad1d` - Row-major `(Q1d * P1d)` matrix expressing derivatives of 447 /// nodal basis functions at quadrature points 448 /// * `qref1d` - Array of length `Q1d` holding the locations of quadrature 449 /// points on the 1D reference element `[-1, 1]` 450 /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on 451 /// the reference element 452 /// 453 /// ``` 454 /// # use libceed::prelude::*; 455 /// # let ceed = libceed::Ceed::default_init(); 456 /// let interp1d = [ 0.62994317, 0.47255875, -0.14950343, 0.04700152, 457 /// -0.07069480, 0.97297619, 0.13253993, -0.03482132, 458 /// -0.03482132, 0.13253993, 0.97297619, -0.07069480, 459 /// 0.04700152, -0.14950343, 0.47255875, 0.62994317]; 460 /// let grad1d = [-2.34183742, 2.78794489, -0.63510411, 0.18899664, 461 /// -0.51670214, -0.48795249, 1.33790510, -0.33325047, 462 // 0.33325047, -1.33790510, 0.48795249, 0.51670214, 463 /// -0.18899664, 0.63510411, -2.78794489, 2.34183742]; 464 /// let qref1d = [-0.86113631, -0.33998104, 0.33998104, 0.86113631]; 465 /// let qweight1d = [ 0.34785485, 0.65214515, 0.65214515, 0.34785485]; 466 /// let b = ceed. 467 /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d).unwrap(); 468 /// ``` 469 pub fn basis_tensor_H1( 470 &self, 471 dim: usize, 472 ncomp: usize, 473 P1d: usize, 474 Q1d: usize, 475 interp1d: &[crate::Scalar], 476 grad1d: &[crate::Scalar], 477 qref1d: &[crate::Scalar], 478 qweight1d: &[crate::Scalar], 479 ) -> Result<Basis> { 480 Basis::create_tensor_H1( 481 self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, 482 ) 483 } 484 485 /// Returns a tensor-product Lagrange basis 486 /// 487 /// # arguments 488 /// 489 /// * `dim` - Topological dimension of element 490 /// * `ncomp` - Number of field components (1 for scalar fields) 491 /// * `P` - Number of Gauss-Lobatto nodes in one dimension. The 492 /// polynomial degree of the resulting `Q_k` element is `k=P-1`. 493 /// * `Q` - Number of quadrature points in one dimension 494 /// * `qmode` - Distribution of the `Q` quadrature points (affects order of 495 /// accuracy for the quadrature) 496 /// 497 /// ``` 498 /// # use libceed::prelude::*; 499 /// # let ceed = libceed::Ceed::default_init(); 500 /// let b = ceed 501 /// .basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss) 502 /// .unwrap(); 503 /// ``` 504 pub fn basis_tensor_H1_Lagrange( 505 &self, 506 dim: usize, 507 ncomp: usize, 508 P: usize, 509 Q: usize, 510 qmode: QuadMode, 511 ) -> Result<Basis> { 512 Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode) 513 } 514 515 /// Returns a tensor-product basis 516 /// 517 /// # arguments 518 /// 519 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 520 /// * `ncomp` - Number of field components (1 for scalar fields) 521 /// * `nnodes` - Total number of nodes 522 /// * `nqpts` - Total number of quadrature points 523 /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of 524 /// nodal basis functions at quadrature points 525 /// * `grad` - Row-major `(nqpts * dim * nnodes)` matrix expressing 526 /// derivatives of nodal basis functions at quadrature points 527 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 528 /// points on the reference element `[-1, 1]` 529 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 530 /// the reference element 531 /// 532 /// ``` 533 /// # use libceed::prelude::*; 534 /// # let ceed = libceed::Ceed::default_init(); 535 /// let interp = [ 536 /// 0.12000000, 537 /// 0.48000000, 538 /// -0.12000000, 539 /// 0.48000000, 540 /// 0.16000000, 541 /// -0.12000000, 542 /// -0.12000000, 543 /// 0.48000000, 544 /// 0.12000000, 545 /// 0.16000000, 546 /// 0.48000000, 547 /// -0.12000000, 548 /// -0.11111111, 549 /// 0.44444444, 550 /// -0.11111111, 551 /// 0.44444444, 552 /// 0.44444444, 553 /// -0.11111111, 554 /// -0.12000000, 555 /// 0.16000000, 556 /// -0.12000000, 557 /// 0.48000000, 558 /// 0.48000000, 559 /// 0.12000000, 560 /// ]; 561 /// let grad = [ 562 /// -1.40000000, 563 /// 1.60000000, 564 /// -0.20000000, 565 /// -0.80000000, 566 /// 0.80000000, 567 /// 0.00000000, 568 /// 0.20000000, 569 /// -1.60000000, 570 /// 1.40000000, 571 /// -0.80000000, 572 /// 0.80000000, 573 /// 0.00000000, 574 /// -0.33333333, 575 /// 0.00000000, 576 /// 0.33333333, 577 /// -1.33333333, 578 /// 1.33333333, 579 /// 0.00000000, 580 /// 0.20000000, 581 /// 0.00000000, 582 /// -0.20000000, 583 /// -2.40000000, 584 /// 2.40000000, 585 /// 0.00000000, 586 /// -1.40000000, 587 /// -0.80000000, 588 /// 0.00000000, 589 /// 1.60000000, 590 /// 0.80000000, 591 /// -0.20000000, 592 /// 0.20000000, 593 /// -2.40000000, 594 /// 0.00000000, 595 /// 0.00000000, 596 /// 2.40000000, 597 /// -0.20000000, 598 /// -0.33333333, 599 /// -1.33333333, 600 /// 0.00000000, 601 /// 0.00000000, 602 /// 1.33333333, 603 /// 0.33333333, 604 /// 0.20000000, 605 /// -0.80000000, 606 /// 0.00000000, 607 /// -1.60000000, 608 /// 0.80000000, 609 /// 1.40000000, 610 /// ]; 611 /// let qref = [ 612 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 613 /// 0.60000000, 614 /// ]; 615 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 616 /// let b = ceed 617 /// .basis_H1( 618 /// ElemTopology::Triangle, 619 /// 1, 620 /// 6, 621 /// 4, 622 /// &interp, 623 /// &grad, 624 /// &qref, 625 /// &qweight, 626 /// ) 627 /// .unwrap(); 628 /// ``` 629 pub fn basis_H1( 630 &self, 631 topo: ElemTopology, 632 ncomp: usize, 633 nnodes: usize, 634 nqpts: usize, 635 interp: &[crate::Scalar], 636 grad: &[crate::Scalar], 637 qref: &[crate::Scalar], 638 qweight: &[crate::Scalar], 639 ) -> Result<Basis> { 640 Basis::create_H1( 641 self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, 642 ) 643 } 644 645 /// Returns a CeedQFunction for evaluating interior (volumetric) terms 646 /// 647 /// # arguments 648 /// 649 /// * `vlength` - Vector length. Caller must ensure that number of 650 /// quadrature points is a multiple of vlength. 651 /// * `f` - Boxed closure to evaluate action at quadrature points. 652 /// 653 /// ``` 654 /// # use libceed::prelude::*; 655 /// # let ceed = libceed::Ceed::default_init(); 656 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 657 /// // Iterate over quadrature points 658 /// v.iter_mut() 659 /// .zip(u.iter().zip(weights.iter())) 660 /// .for_each(|(v, (u, w))| *v = u * w); 661 /// 662 /// // Return clean error code 663 /// 0 664 /// }; 665 /// 666 /// let qf = ceed.q_function_interior(1, Box::new(user_f)).unwrap(); 667 /// ``` 668 pub fn q_function_interior( 669 &self, 670 vlength: usize, 671 f: Box<qfunction::QFunctionUserClosure>, 672 ) -> Result<QFunction> { 673 QFunction::create(self, vlength, f) 674 } 675 676 /// Returns a CeedQFunction for evaluating interior (volumetric) terms 677 /// created by name 678 /// 679 /// ``` 680 /// # use libceed::prelude::*; 681 /// # let ceed = libceed::Ceed::default_init(); 682 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 683 /// ``` 684 pub fn q_function_interior_by_name(&self, name: &str) -> Result<QFunctionByName> { 685 QFunctionByName::create(self, name) 686 } 687 688 /// Returns a Operator and associate a QFunction. A Basis and 689 /// ElemRestriction can be associated with QFunction fields with 690 /// set_field(). 691 /// 692 /// * `qf` - QFunction defining the action of the operator at quadrature 693 /// points 694 /// * `dqf` - QFunction defining the action of the Jacobian of the qf (or 695 /// qfunction_none) 696 /// * `dqfT` - QFunction defining the action of the transpose of the 697 /// Jacobian of the qf (or qfunction_none) 698 /// 699 /// ``` 700 /// # use libceed::prelude::*; 701 /// # let ceed = libceed::Ceed::default_init(); 702 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild").unwrap(); 703 /// let op = ceed 704 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None) 705 /// .unwrap(); 706 /// ``` 707 pub fn operator<'b>( 708 &self, 709 qf: impl Into<QFunctionOpt<'b>>, 710 dqf: impl Into<QFunctionOpt<'b>>, 711 dqfT: impl Into<QFunctionOpt<'b>>, 712 ) -> Result<Operator> { 713 Operator::create(self, qf, dqf, dqfT) 714 } 715 716 /// Returns an Operator that composes the action of several Operators 717 /// 718 /// ``` 719 /// # use libceed::prelude::*; 720 /// # let ceed = libceed::Ceed::default_init(); 721 /// let op = ceed.composite_operator().unwrap(); 722 /// ``` 723 pub fn composite_operator(&self) -> Result<CompositeOperator> { 724 CompositeOperator::create(self) 725 } 726 } 727 728 // ----------------------------------------------------------------------------- 729 // Tests 730 // ----------------------------------------------------------------------------- 731 #[cfg(test)] 732 mod tests { 733 use super::*; 734 735 fn ceed_t501() -> Result<i32> { 736 let resource = "/cpu/self/ref/blocked"; 737 let ceed = Ceed::init(resource); 738 let nelem = 4; 739 let p = 3; 740 let q = 4; 741 let ndofs = p * nelem - nelem + 1; 742 743 // Vectors 744 let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 745 let mut qdata = ceed.vector(nelem * q)?; 746 qdata.set_value(0.0)?; 747 let mut u = ceed.vector(ndofs)?; 748 u.set_value(1.0)?; 749 let mut v = ceed.vector(ndofs)?; 750 v.set_value(0.0)?; 751 752 // Restrictions 753 let mut indx: Vec<i32> = vec![0; 2 * nelem]; 754 for i in 0..nelem { 755 indx[2 * i + 0] = i as i32; 756 indx[2 * i + 1] = (i + 1) as i32; 757 } 758 let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?; 759 let mut indu: Vec<i32> = vec![0; p * nelem]; 760 for i in 0..nelem { 761 indu[p * i + 0] = i as i32; 762 indu[p * i + 1] = (i + 1) as i32; 763 indu[p * i + 2] = (i + 2) as i32; 764 } 765 let ru = ceed.elem_restriction(nelem, 3, 1, 1, ndofs, MemType::Host, &indu)?; 766 let strides: [i32; 3] = [1, q as i32, q as i32]; 767 let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?; 768 769 // Bases 770 let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 771 let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 772 773 // Build quadrature data 774 let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 775 ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 776 .field("dx", &rx, &bx, VectorOpt::Active)? 777 .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 778 .field("qdata", &rq, BasisOpt::Collocated, VectorOpt::Active)? 779 .apply(&x, &mut qdata)?; 780 781 // Mass operator 782 let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 783 let op_mass = ceed 784 .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 785 .field("u", &ru, &bu, VectorOpt::Active)? 786 .field("qdata", &rq, BasisOpt::Collocated, &qdata)? 787 .field("v", &ru, &bu, VectorOpt::Active)?; 788 789 v.set_value(0.0)?; 790 op_mass.apply(&u, &mut v)?; 791 792 // Check 793 let sum: Scalar = v.view().iter().sum(); 794 assert!( 795 (sum - 2.0).abs() < 1e-15, 796 "Incorrect interval length computed" 797 ); 798 Ok(0) 799 } 800 801 #[test] 802 fn test_ceed_t501() { 803 assert!(ceed_t501().is_ok()); 804 } 805 } 806 807 // ----------------------------------------------------------------------------- 808