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