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