1 // Copyright (c) 2017-2024, 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: usize = 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<F>(ceed_ptr: F, ierr: i32) -> Result<i32> 165 where 166 F: FnOnce() -> bind_ceed::Ceed, 167 { 168 // Return early if code is clean 169 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 170 return Ok(ierr); 171 } 172 // Retrieve error message 173 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 174 let c_str = unsafe { 175 bind_ceed::CeedGetErrorMessage(ceed_ptr(), &mut ptr); 176 std::ffi::CStr::from_ptr(ptr) 177 }; 178 let message = c_str.to_string_lossy().to_string(); 179 // Panic if negative code, otherwise return error 180 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 181 panic!("{}", message); 182 } 183 Err(Error { message }) 184 } 185 186 // ----------------------------------------------------------------------------- 187 // Ceed error handler 188 // ----------------------------------------------------------------------------- 189 pub enum ErrorHandler { 190 ErrorAbort, 191 ErrorExit, 192 ErrorReturn, 193 ErrorStore, 194 } 195 196 // ----------------------------------------------------------------------------- 197 // Ceed context wrapper 198 // ----------------------------------------------------------------------------- 199 /// A Ceed is a library context representing control of a logical hardware 200 /// resource. 201 #[derive(Debug)] 202 pub struct Ceed { 203 ptr: bind_ceed::Ceed, 204 } 205 206 // ----------------------------------------------------------------------------- 207 // Destructor 208 // ----------------------------------------------------------------------------- 209 impl Drop for Ceed { 210 fn drop(&mut self) { 211 unsafe { 212 bind_ceed::CeedDestroy(&mut self.ptr); 213 } 214 } 215 } 216 217 // ----------------------------------------------------------------------------- 218 // Cloning 219 // ----------------------------------------------------------------------------- 220 impl Clone for Ceed { 221 /// Perform a shallow clone of a Ceed context 222 /// 223 /// ``` 224 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 225 /// let ceed_clone = ceed.clone(); 226 /// 227 /// println!(" original:{} \n clone:{}", ceed, ceed_clone); 228 /// ``` 229 fn clone(&self) -> Self { 230 let mut ptr_clone = std::ptr::null_mut(); 231 self.check_error(unsafe { bind_ceed::CeedReferenceCopy(self.ptr, &mut ptr_clone) }) 232 .expect("failed to clone Ceed"); 233 Self { ptr: ptr_clone } 234 } 235 } 236 237 // ----------------------------------------------------------------------------- 238 // Display 239 // ----------------------------------------------------------------------------- 240 impl fmt::Display for Ceed { 241 /// View a Ceed 242 /// 243 /// ``` 244 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 245 /// println!("{}", ceed); 246 /// ``` 247 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 248 let mut ptr = std::ptr::null_mut(); 249 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 250 let cstring = unsafe { 251 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 252 bind_ceed::CeedView(self.ptr, file); 253 bind_ceed::fclose(file); 254 CString::from_raw(ptr) 255 }; 256 cstring.to_string_lossy().fmt(f) 257 } 258 } 259 260 static REGISTER: Once = Once::new(); 261 262 // ----------------------------------------------------------------------------- 263 // Object constructors 264 // ----------------------------------------------------------------------------- 265 impl Ceed { 266 #[cfg_attr(feature = "katexit", katexit::katexit)] 267 /// Returns a Ceed context initialized with the specified resource 268 /// 269 /// # arguments 270 /// 271 /// * `resource` - Resource to use, e.g., "/cpu/self" 272 /// 273 /// ``` 274 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 275 /// ``` 276 pub fn init(resource: &str) -> Self { 277 Ceed::init_with_error_handler(resource, ErrorHandler::ErrorStore) 278 } 279 280 /// Returns a Ceed context initialized with the specified resource 281 /// 282 /// # arguments 283 /// 284 /// * `resource` - Resource to use, e.g., "/cpu/self" 285 /// 286 /// ``` 287 /// let ceed = libceed::Ceed::init_with_error_handler( 288 /// "/cpu/self/ref/serial", 289 /// libceed::ErrorHandler::ErrorAbort, 290 /// ); 291 /// ``` 292 pub fn init_with_error_handler(resource: &str, handler: ErrorHandler) -> Self { 293 REGISTER.call_once(|| unsafe { 294 bind_ceed::CeedRegisterAll(); 295 bind_ceed::CeedQFunctionRegisterAll(); 296 }); 297 298 // Convert to C string 299 let c_resource = CString::new(resource).expect("CString::new failed"); 300 301 // Get error handler pointer 302 let eh = match handler { 303 ErrorHandler::ErrorAbort => bind_ceed::CeedErrorAbort, 304 ErrorHandler::ErrorExit => bind_ceed::CeedErrorExit, 305 ErrorHandler::ErrorReturn => bind_ceed::CeedErrorReturn, 306 ErrorHandler::ErrorStore => bind_ceed::CeedErrorStore, 307 }; 308 309 // Call to libCEED 310 let mut ptr = std::ptr::null_mut(); 311 let mut ierr = unsafe { bind_ceed::CeedInit(c_resource.as_ptr() as *const i8, &mut ptr) }; 312 if ierr != 0 { 313 panic!("Error initializing backend resource: {}", resource) 314 } 315 ierr = unsafe { bind_ceed::CeedSetErrorHandler(ptr, Some(eh)) }; 316 let ceed = Ceed { ptr }; 317 ceed.check_error(ierr).unwrap(); 318 ceed 319 } 320 321 /// Default initializer for testing 322 #[doc(hidden)] 323 pub fn default_init() -> Self { 324 // Convert to C string 325 let resource = "/cpu/self/ref/serial"; 326 crate::Ceed::init(resource) 327 } 328 329 /// Internal error checker 330 #[doc(hidden)] 331 fn check_error(&self, ierr: i32) -> Result<i32> { 332 // Return early if code is clean 333 if ierr == bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 334 return Ok(ierr); 335 } 336 // Retrieve error message 337 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 338 let c_str = unsafe { 339 bind_ceed::CeedGetErrorMessage(self.ptr, &mut ptr); 340 std::ffi::CStr::from_ptr(ptr) 341 }; 342 let message = c_str.to_string_lossy().to_string(); 343 // Panic if negative code, otherwise return error 344 if ierr < bind_ceed::CeedErrorType_CEED_ERROR_SUCCESS { 345 panic!("{}", message); 346 } 347 Err(Error { message }) 348 } 349 350 /// Returns full resource name for a Ceed context 351 /// 352 /// ``` 353 /// let ceed = libceed::Ceed::init("/cpu/self/ref/serial"); 354 /// let resource = ceed.resource(); 355 /// 356 /// assert_eq!(resource, "/cpu/self/ref/serial".to_string()) 357 /// ``` 358 pub fn resource(&self) -> String { 359 let mut ptr: *const std::os::raw::c_char = std::ptr::null_mut(); 360 let c_str = unsafe { 361 bind_ceed::CeedGetResource(self.ptr, &mut ptr); 362 std::ffi::CStr::from_ptr(ptr) 363 }; 364 c_str.to_string_lossy().to_string() 365 } 366 367 /// Returns a Vector of the specified length (does not allocate memory) 368 /// 369 /// # arguments 370 /// 371 /// * `n` - Length of vector 372 /// 373 /// ``` 374 /// # use libceed::prelude::*; 375 /// # fn main() -> libceed::Result<()> { 376 /// # let ceed = libceed::Ceed::default_init(); 377 /// let vec = ceed.vector(10)?; 378 /// # Ok(()) 379 /// # } 380 /// ``` 381 pub fn vector<'a>(&self, n: usize) -> Result<Vector<'a>> { 382 Vector::create(self, n) 383 } 384 385 /// Create a Vector initialized with the data (copied) from a slice 386 /// 387 /// # arguments 388 /// 389 /// * `slice` - Slice containing data 390 /// 391 /// ``` 392 /// # use libceed::prelude::*; 393 /// # fn main() -> libceed::Result<()> { 394 /// # let ceed = libceed::Ceed::default_init(); 395 /// let vec = ceed.vector_from_slice(&[1., 2., 3.])?; 396 /// assert_eq!(vec.length(), 3); 397 /// # Ok(()) 398 /// # } 399 /// ``` 400 pub fn vector_from_slice<'a>(&self, slice: &[crate::Scalar]) -> Result<Vector<'a>> { 401 Vector::from_slice(self, slice) 402 } 403 404 /// Returns an ElemRestriction, $\mathcal{E}$, which extracts the degrees of 405 /// freedom for each element from the local vector into the element vector 406 /// or assembles contributions from each element in the element vector to 407 /// the local vector 408 /// 409 /// # arguments 410 /// 411 /// * `nelem` - Number of elements described in the offsets array 412 /// * `elemsize` - Size (number of "nodes") per element 413 /// * `ncomp` - Number of field components per interpolation node (1 414 /// for scalar fields) 415 /// * `compstride` - Stride between components for the same Lvector "node". 416 /// Data for node `i`, component `j`, element `k` can be 417 /// found in the Lvector at index 418 /// `offsets[i + k*elemsize] + j*compstride`. 419 /// * `lsize` - The size of the Lvector. This vector may be larger 420 /// than the elements and fields given by this 421 /// restriction. 422 /// * `mtype` - Memory type of the offsets array, see CeedMemType 423 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 424 /// ordered list of the offsets (into the input CeedVector) 425 /// for the unknowns corresponding to element `i`, where 426 /// `0 <= i < nelem`. All offsets must be in the range 427 /// `[0, lsize - 1]`. 428 /// 429 /// ``` 430 /// # use libceed::prelude::*; 431 /// # fn main() -> libceed::Result<()> { 432 /// # let ceed = libceed::Ceed::default_init(); 433 /// let nelem = 3; 434 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 435 /// for i in 0..nelem { 436 /// ind[2 * i + 0] = i as i32; 437 /// ind[2 * i + 1] = (i + 1) as i32; 438 /// } 439 /// let r = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind)?; 440 /// # Ok(()) 441 /// # } 442 /// ``` 443 pub fn elem_restriction<'a>( 444 &self, 445 nelem: usize, 446 elemsize: usize, 447 ncomp: usize, 448 compstride: usize, 449 lsize: usize, 450 mtype: MemType, 451 offsets: &[i32], 452 ) -> Result<ElemRestriction<'a>> { 453 ElemRestriction::create( 454 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, 455 ) 456 } 457 458 /// Returns an oriented ElemRestriction, $\mathcal{E}$, which extracts the 459 /// degrees of freedom for each element from the local vector into the 460 /// element vector or assembles contributions from each element in the 461 /// element vector to the local vector 462 /// 463 /// # arguments 464 /// 465 /// * `nelem` - Number of elements described in the offsets array 466 /// * `elemsize` - Size (number of "nodes") per element 467 /// * `ncomp` - Number of field components per interpolation node (1 468 /// for scalar fields) 469 /// * `compstride` - Stride between components for the same Lvector "node". 470 /// Data for node `i`, component `j`, element `k` can be 471 /// found in the Lvector at index 472 /// `offsets[i + k*elemsize] + j*compstride`. 473 /// * `lsize` - The size of the Lvector. This vector may be larger 474 /// than the elements and fields given by this 475 /// restriction. 476 /// * `mtype` - Memory type of the offsets array, see CeedMemType 477 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 478 /// ordered list of the offsets (into the input CeedVector) 479 /// for the unknowns corresponding to element `i`, where 480 /// `0 <= i < nelem`. All offsets must be in the range 481 /// `[0, lsize - 1]`. 482 /// * `orients` - Array of shape `[nelem, elemsize]`. Row `i` holds the 483 /// ordered list of the orientations for the unknowns 484 /// corresponding to element `i`, with bool `false` used 485 /// for positively oriented and `true` to flip the 486 /// orientation. 487 /// 488 /// ``` 489 /// # use libceed::prelude::*; 490 /// # fn main() -> libceed::Result<()> { 491 /// # let ceed = libceed::Ceed::default_init(); 492 /// let nelem = 3; 493 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 494 /// let mut orients: Vec<bool> = vec![false; 2 * nelem]; 495 /// for i in 0..nelem { 496 /// ind[2 * i + 0] = i as i32; 497 /// ind[2 * i + 1] = (i + 1) as i32; 498 /// orients[2 * i + 0] = (i % 2) > 0; // flip the dofs on element 1, 3, ... 499 /// orients[2 * i + 1] = (i % 2) > 0; 500 /// } 501 /// let r = 502 /// ceed.oriented_elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &ind, &orients)?; 503 /// # Ok(()) 504 /// # } 505 /// ``` 506 pub fn oriented_elem_restriction<'a>( 507 &self, 508 nelem: usize, 509 elemsize: usize, 510 ncomp: usize, 511 compstride: usize, 512 lsize: usize, 513 mtype: MemType, 514 offsets: &[i32], 515 orients: &[bool], 516 ) -> Result<ElemRestriction<'a>> { 517 ElemRestriction::create_oriented( 518 self, nelem, elemsize, ncomp, compstride, lsize, mtype, offsets, orients, 519 ) 520 } 521 522 /// Returns a curl-oriented ElemRestriction, $\mathcal{E}$, which extracts 523 /// the degrees of freedom for each element from the local vector into 524 /// the element vector or assembles contributions from each element in 525 /// the element vector to the local vector 526 /// 527 /// # arguments 528 /// 529 /// * `nelem` - Number of elements described in the offsets array 530 /// * `elemsize` - Size (number of "nodes") per element 531 /// * `ncomp` - Number of field components per interpolation node (1 532 /// for scalar fields) 533 /// * `compstride` - Stride between components for the same Lvector "node". 534 /// Data for node `i`, component `j`, element `k` can be 535 /// found in the Lvector at index 536 /// `offsets[i + k*elemsize] + j*compstride`. 537 /// * `lsize` - The size of the Lvector. This vector may be larger 538 /// than the elements and fields given by this 539 /// restriction. 540 /// * `mtype` - Memory type of the offsets array, see CeedMemType 541 /// * `offsets` - Array of shape `[nelem, elemsize]`. Row `i` holds the 542 /// ordered list of the offsets (into the input CeedVector) 543 /// for the unknowns corresponding to element `i`, where 544 /// `0 <= i < nelem`. All offsets must be in the range 545 /// `[0, lsize - 1]`. 546 /// * `curlorients` - Array of shape `[nelem, 3 * elemsize]`. Row `i` holds 547 /// a row-major tridiagonal matrix (`curlorients[i, 0] = 548 /// curlorients[i, 3 * elemsize - 1] = 0`, where 549 /// `0 <= i < nelem`) which is applied to the element 550 /// unknowns upon restriction. 551 /// 552 /// ``` 553 /// # use libceed::prelude::*; 554 /// # fn main() -> libceed::Result<()> { 555 /// # let ceed = libceed::Ceed::default_init(); 556 /// let nelem = 3; 557 /// let mut ind: Vec<i32> = vec![0; 2 * nelem]; 558 /// let mut curlorients: Vec<i8> = vec![0; 3 * 2 * nelem]; 559 /// for i in 0..nelem { 560 /// ind[2 * i + 0] = i as i32; 561 /// ind[2 * i + 1] = (i + 1) as i32; 562 /// curlorients[3 * 2 * i] = 0 as i8; 563 /// curlorients[3 * 2 * (i + 1) - 1] = 0 as i8; 564 /// if (i % 2 > 0) { 565 /// // T = [0 -1] 566 /// // [-1 0] 567 /// curlorients[3 * 2 * i + 1] = 0 as i8; 568 /// curlorients[3 * 2 * i + 2] = -1 as i8; 569 /// curlorients[3 * 2 * i + 3] = -1 as i8; 570 /// curlorients[3 * 2 * i + 4] = 0 as i8; 571 /// } else { 572 /// // T = I 573 /// curlorients[3 * 2 * i + 1] = 1 as i8; 574 /// curlorients[3 * 2 * i + 2] = 0 as i8; 575 /// curlorients[3 * 2 * i + 3] = 0 as i8; 576 /// curlorients[3 * 2 * i + 4] = 1 as i8; 577 /// } 578 /// } 579 /// let r = ceed.curl_oriented_elem_restriction( 580 /// nelem, 581 /// 2, 582 /// 1, 583 /// 1, 584 /// nelem + 1, 585 /// MemType::Host, 586 /// &ind, 587 /// &curlorients, 588 /// )?; 589 /// # Ok(()) 590 /// # } 591 /// ``` 592 pub fn curl_oriented_elem_restriction<'a>( 593 &self, 594 nelem: usize, 595 elemsize: usize, 596 ncomp: usize, 597 compstride: usize, 598 lsize: usize, 599 mtype: MemType, 600 offsets: &[i32], 601 curlorients: &[i8], 602 ) -> Result<ElemRestriction<'a>> { 603 ElemRestriction::create_curl_oriented( 604 self, 605 nelem, 606 elemsize, 607 ncomp, 608 compstride, 609 lsize, 610 mtype, 611 offsets, 612 curlorients, 613 ) 614 } 615 616 /// Returns an ElemRestriction, $\mathcal{E}$, from an local vector to 617 /// an element vector where data can be indexed from the `strides` array 618 /// 619 /// # arguments 620 /// 621 /// * `nelem` - Number of elements described in the offsets array 622 /// * `elemsize` - Size (number of "nodes") per element 623 /// * `ncomp` - Number of field components per interpolation node (1 624 /// for scalar fields) 625 /// * `lsize` - The size of the Lvector. This vector may be larger 626 /// than the elements and fields given by this restriction. 627 /// * `strides` - Array for strides between `[nodes, components, elements]`. 628 /// Data for node `i`, component `j`, element `k` can be 629 /// found in the Lvector at index 630 /// `i*strides[0] + j*strides[1] + k*strides[2]`. 631 /// CEED_STRIDES_BACKEND may be used with vectors created 632 /// by a Ceed backend. 633 /// 634 /// ``` 635 /// # use libceed::prelude::*; 636 /// # fn main() -> libceed::Result<()> { 637 /// # let ceed = libceed::Ceed::default_init(); 638 /// let nelem = 3; 639 /// let strides: [i32; 3] = [1, 2, 2]; 640 /// let r = ceed.strided_elem_restriction(nelem, 2, 1, nelem * 2, strides)?; 641 /// # Ok(()) 642 /// # } 643 /// ``` 644 pub fn strided_elem_restriction<'a>( 645 &self, 646 nelem: usize, 647 elemsize: usize, 648 ncomp: usize, 649 lsize: usize, 650 strides: [i32; 3], 651 ) -> Result<ElemRestriction<'a>> { 652 ElemRestriction::create_strided(self, nelem, elemsize, ncomp, lsize, strides) 653 } 654 655 /// Returns an $H^1$ tensor-product Basis 656 /// 657 /// # arguments 658 /// 659 /// * `dim` - Topological dimension of element 660 /// * `ncomp` - Number of field components (1 for scalar fields) 661 /// * `P1d` - Number of Gauss-Lobatto nodes in one dimension. The 662 /// polynomial degree of the resulting `Q_k` element is 663 /// `k=P-1`. 664 /// * `Q1d` - Number of quadrature points in one dimension 665 /// * `interp1d` - Row-major `(Q1d * P1d)` matrix expressing the values of 666 /// nodal basis functions at quadrature points 667 /// * `grad1d` - Row-major `(Q1d * P1d)` matrix expressing derivatives of 668 /// nodal basis functions at quadrature points 669 /// * `qref1d` - Array of length `Q1d` holding the locations of quadrature 670 /// points on the 1D reference element `[-1, 1]` 671 /// * `qweight1d` - Array of length `Q1d` holding the quadrature weights on 672 /// the reference element 673 /// 674 /// ``` 675 /// # use libceed::prelude::*; 676 /// # fn main() -> libceed::Result<()> { 677 /// # let ceed = libceed::Ceed::default_init(); 678 /// let interp1d = [ 0.62994317, 0.47255875, -0.14950343, 0.04700152, 679 /// -0.07069480, 0.97297619, 0.13253993, -0.03482132, 680 /// -0.03482132, 0.13253993, 0.97297619, -0.07069480, 681 /// 0.04700152, -0.14950343, 0.47255875, 0.62994317]; 682 /// let grad1d = [-2.34183742, 2.78794489, -0.63510411, 0.18899664, 683 /// -0.51670214, -0.48795249, 1.33790510, -0.33325047, 684 // 0.33325047, -1.33790510, 0.48795249, 0.51670214, 685 /// -0.18899664, 0.63510411, -2.78794489, 2.34183742]; 686 /// let qref1d = [-0.86113631, -0.33998104, 0.33998104, 0.86113631]; 687 /// let qweight1d = [ 0.34785485, 0.65214515, 0.65214515, 0.34785485]; 688 /// let b = ceed. 689 /// basis_tensor_H1(2, 1, 4, 4, &interp1d, &grad1d, &qref1d, &qweight1d)?; 690 /// # Ok(()) 691 /// # } 692 /// ``` 693 pub fn basis_tensor_H1<'a>( 694 &self, 695 dim: usize, 696 ncomp: usize, 697 P1d: usize, 698 Q1d: usize, 699 interp1d: &[crate::Scalar], 700 grad1d: &[crate::Scalar], 701 qref1d: &[crate::Scalar], 702 qweight1d: &[crate::Scalar], 703 ) -> Result<Basis<'a>> { 704 Basis::create_tensor_H1( 705 self, dim, ncomp, P1d, Q1d, interp1d, grad1d, qref1d, qweight1d, 706 ) 707 } 708 709 /// Returns an $H^1$ Lagrange tensor-product Basis 710 /// 711 /// # arguments 712 /// 713 /// * `dim` - Topological dimension of element 714 /// * `ncomp` - Number of field components (1 for scalar fields) 715 /// * `P` - Number of Gauss-Lobatto nodes in one dimension. The 716 /// polynomial degree of the resulting `Q_k` element is `k=P-1`. 717 /// * `Q` - Number of quadrature points in one dimension 718 /// * `qmode` - Distribution of the `Q` quadrature points (affects order of 719 /// accuracy for the quadrature) 720 /// 721 /// ``` 722 /// # use libceed::prelude::*; 723 /// # fn main() -> libceed::Result<()> { 724 /// # let ceed = libceed::Ceed::default_init(); 725 /// let b = ceed.basis_tensor_H1_Lagrange(2, 1, 3, 4, QuadMode::Gauss)?; 726 /// # Ok(()) 727 /// # } 728 /// ``` 729 pub fn basis_tensor_H1_Lagrange<'a>( 730 &self, 731 dim: usize, 732 ncomp: usize, 733 P: usize, 734 Q: usize, 735 qmode: QuadMode, 736 ) -> Result<Basis<'a>> { 737 Basis::create_tensor_H1_Lagrange(self, dim, ncomp, P, Q, qmode) 738 } 739 740 /// Returns an $H-1$ Basis 741 /// 742 /// # arguments 743 /// 744 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 745 /// * `ncomp` - Number of field components (1 for scalar fields) 746 /// * `nnodes` - Total number of nodes 747 /// * `nqpts` - Total number of quadrature points 748 /// * `interp` - Row-major `(nqpts * nnodes)` matrix expressing the values of 749 /// nodal basis functions at quadrature points 750 /// * `grad` - Row-major `(dim * nqpts * nnodes)` matrix expressing 751 /// derivatives of nodal basis functions at quadrature points 752 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 753 /// points on the reference element `[-1, 1]` 754 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 755 /// the reference element 756 /// 757 /// ``` 758 /// # use libceed::prelude::*; 759 /// # fn main() -> libceed::Result<()> { 760 /// # let ceed = libceed::Ceed::default_init(); 761 /// let interp = [ 762 /// 0.12000000, 763 /// 0.48000000, 764 /// -0.12000000, 765 /// 0.48000000, 766 /// 0.16000000, 767 /// -0.12000000, 768 /// -0.12000000, 769 /// 0.48000000, 770 /// 0.12000000, 771 /// 0.16000000, 772 /// 0.48000000, 773 /// -0.12000000, 774 /// -0.11111111, 775 /// 0.44444444, 776 /// -0.11111111, 777 /// 0.44444444, 778 /// 0.44444444, 779 /// -0.11111111, 780 /// -0.12000000, 781 /// 0.16000000, 782 /// -0.12000000, 783 /// 0.48000000, 784 /// 0.48000000, 785 /// 0.12000000, 786 /// ]; 787 /// let grad = [ 788 /// -1.40000000, 789 /// 1.60000000, 790 /// -0.20000000, 791 /// -0.80000000, 792 /// 0.80000000, 793 /// 0.00000000, 794 /// 0.20000000, 795 /// -1.60000000, 796 /// 1.40000000, 797 /// -0.80000000, 798 /// 0.80000000, 799 /// 0.00000000, 800 /// -0.33333333, 801 /// 0.00000000, 802 /// 0.33333333, 803 /// -1.33333333, 804 /// 1.33333333, 805 /// 0.00000000, 806 /// 0.20000000, 807 /// 0.00000000, 808 /// -0.20000000, 809 /// -2.40000000, 810 /// 2.40000000, 811 /// 0.00000000, 812 /// -1.40000000, 813 /// -0.80000000, 814 /// 0.00000000, 815 /// 1.60000000, 816 /// 0.80000000, 817 /// -0.20000000, 818 /// 0.20000000, 819 /// -2.40000000, 820 /// 0.00000000, 821 /// 0.00000000, 822 /// 2.40000000, 823 /// -0.20000000, 824 /// -0.33333333, 825 /// -1.33333333, 826 /// 0.00000000, 827 /// 0.00000000, 828 /// 1.33333333, 829 /// 0.33333333, 830 /// 0.20000000, 831 /// -0.80000000, 832 /// 0.00000000, 833 /// -1.60000000, 834 /// 0.80000000, 835 /// 1.40000000, 836 /// ]; 837 /// let qref = [ 838 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 839 /// 0.60000000, 840 /// ]; 841 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 842 /// let b = ceed.basis_H1( 843 /// ElemTopology::Triangle, 844 /// 1, 845 /// 6, 846 /// 4, 847 /// &interp, 848 /// &grad, 849 /// &qref, 850 /// &qweight, 851 /// )?; 852 /// # Ok(()) 853 /// # } 854 /// ``` 855 pub fn basis_H1<'a>( 856 &self, 857 topo: ElemTopology, 858 ncomp: usize, 859 nnodes: usize, 860 nqpts: usize, 861 interp: &[crate::Scalar], 862 grad: &[crate::Scalar], 863 qref: &[crate::Scalar], 864 qweight: &[crate::Scalar], 865 ) -> Result<Basis<'a>> { 866 Basis::create_H1( 867 self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight, 868 ) 869 } 870 871 /// Returns an $H(div)$ Basis 872 /// 873 /// # arguments 874 /// 875 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 876 /// * `ncomp` - Number of field components (1 for scalar fields) 877 /// * `nnodes` - Total number of nodes 878 /// * `nqpts` - Total number of quadrature points 879 /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the 880 /// values of basis functions at quadrature points 881 /// * `div` - Row-major `(nqpts * nnodes)` matrix expressing the 882 /// divergence of basis functions at quadrature points 883 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 884 /// points on the reference element `[-1, 1]` 885 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 886 /// the reference element 887 /// 888 /// ``` 889 /// # use libceed::prelude::*; 890 /// # fn main() -> libceed::Result<()> { 891 /// # let ceed = libceed::Ceed::default_init(); 892 /// let interp = [ 893 /// 0.00000000, 894 /// -1.57735027, 895 /// 0.57735027, 896 /// 0.00000000, 897 /// 0.00000000, 898 /// -1.57735027, 899 /// 0.57735027, 900 /// 0.00000000, 901 /// 0.00000000, 902 /// -1.57735027, 903 /// 0.57735027, 904 /// 0.00000000, 905 /// 0.00000000, 906 /// -0.42264973, 907 /// -0.57735027, 908 /// 0.00000000, 909 /// 0.42264973, 910 /// 0.00000000, 911 /// 0.00000000, 912 /// 0.57735027, 913 /// 0.42264973, 914 /// 0.00000000, 915 /// 0.00000000, 916 /// 0.57735027, 917 /// 1.57735027, 918 /// 0.00000000, 919 /// 0.00000000, 920 /// -0.57735027, 921 /// 1.57735027, 922 /// 0.00000000, 923 /// 0.00000000, 924 /// -0.57735027, 925 /// ]; 926 /// let div = [ 927 /// -1.00000000, 928 /// 1.00000000, 929 /// -1.00000000, 930 /// 1.00000000, 931 /// -1.00000000, 932 /// 1.00000000, 933 /// -1.00000000, 934 /// 1.00000000, 935 /// -1.00000000, 936 /// 1.00000000, 937 /// -1.00000000, 938 /// 1.00000000, 939 /// -1.00000000, 940 /// 1.00000000, 941 /// -1.00000000, 942 /// 1.00000000, 943 /// ]; 944 /// let qref = [ 945 /// 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 0.57735026, 946 /// 0.57735026, 947 /// ]; 948 /// let qweight = [1.00000000, 1.00000000, 1.00000000, 1.00000000]; 949 /// let b = ceed.basis_Hdiv(ElemTopology::Quad, 1, 4, 4, &interp, &div, &qref, &qweight)?; 950 /// # Ok(()) 951 /// # } 952 /// ``` 953 pub fn basis_Hdiv<'a>( 954 &self, 955 topo: ElemTopology, 956 ncomp: usize, 957 nnodes: usize, 958 nqpts: usize, 959 interp: &[crate::Scalar], 960 div: &[crate::Scalar], 961 qref: &[crate::Scalar], 962 qweight: &[crate::Scalar], 963 ) -> Result<Basis<'a>> { 964 Basis::create_Hdiv(self, topo, ncomp, nnodes, nqpts, interp, div, qref, qweight) 965 } 966 967 /// Returns an $H(curl)$ Basis 968 /// 969 /// # arguments 970 /// 971 /// * `topo` - Topology of element, e.g. hypercube, simplex, ect 972 /// * `ncomp` - Number of field components (1 for scalar fields) 973 /// * `nnodes` - Total number of nodes 974 /// * `nqpts` - Total number of quadrature points 975 /// * `interp` - Row-major `(dim * nqpts * nnodes)` matrix expressing the 976 /// values of basis functions at quadrature points 977 /// * `curl` - Row-major `(curl_comp * nqpts * nnodes)`, `curl_comp = 1 if 978 /// dim < 3 else dim` matrix expressing the curl of basis 979 /// functions at quadrature points 980 /// * `qref` - Array of length `nqpts` holding the locations of quadrature 981 /// points on the reference element `[-1, 1]` 982 /// * `qweight` - Array of length `nqpts` holding the quadrature weights on 983 /// the reference element 984 /// 985 /// ``` 986 /// # use libceed::prelude::*; 987 /// # fn main() -> libceed::Result<()> { 988 /// # let ceed = libceed::Ceed::default_init(); 989 /// let interp = [ 990 /// -0.20000000, 991 /// 0.20000000, 992 /// 0.80000000, 993 /// -0.20000000, 994 /// 0.20000000, 995 /// 0.80000000, 996 /// -0.33333333, 997 /// 0.33333333, 998 /// 0.66666667, 999 /// -0.60000000, 1000 /// 0.60000000, 1001 /// 0.40000000, 1002 /// 0.20000000, 1003 /// 0.80000000, 1004 /// 0.20000000, 1005 /// 0.60000000, 1006 /// 0.40000000, 1007 /// 0.60000000, 1008 /// 0.33333333, 1009 /// 0.66666667, 1010 /// 0.33333333, 1011 /// 0.20000000, 1012 /// 0.80000000, 1013 /// 0.20000000, 1014 /// ]; 1015 /// let curl = [ 1016 /// 2.00000000, 1017 /// -2.00000000, 1018 /// -2.00000000, 1019 /// 2.00000000, 1020 /// -2.00000000, 1021 /// -2.00000000, 1022 /// 2.00000000, 1023 /// -2.00000000, 1024 /// -2.00000000, 1025 /// 2.00000000, 1026 /// -2.00000000, 1027 /// -2.00000000, 1028 /// ]; 1029 /// let qref = [ 1030 /// 0.20000000, 0.60000000, 0.33333333, 0.20000000, 0.20000000, 0.20000000, 0.33333333, 1031 /// 0.60000000, 1032 /// ]; 1033 /// let qweight = [0.26041667, 0.26041667, -0.28125000, 0.26041667]; 1034 /// let b = ceed.basis_Hcurl( 1035 /// ElemTopology::Triangle, 1036 /// 1, 1037 /// 3, 1038 /// 4, 1039 /// &interp, 1040 /// &curl, 1041 /// &qref, 1042 /// &qweight, 1043 /// )?; 1044 /// # Ok(()) 1045 /// # } 1046 /// ``` 1047 pub fn basis_Hcurl<'a>( 1048 &self, 1049 topo: ElemTopology, 1050 ncomp: usize, 1051 nnodes: usize, 1052 nqpts: usize, 1053 interp: &[crate::Scalar], 1054 curl: &[crate::Scalar], 1055 qref: &[crate::Scalar], 1056 qweight: &[crate::Scalar], 1057 ) -> Result<Basis<'a>> { 1058 Basis::create_Hcurl( 1059 self, topo, ncomp, nnodes, nqpts, interp, curl, qref, qweight, 1060 ) 1061 } 1062 1063 /// Returns a QFunction for evaluating interior (volumetric) terms 1064 /// of a weak form corresponding to the $L^2$ inner product 1065 /// 1066 /// $$ 1067 /// \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), 1068 /// $$ 1069 /// 1070 /// where $v \cdot f_0$ represents contraction over fields and $\nabla v : f_1$ 1071 /// represents contraction over both fields and spatial dimensions. 1072 /// 1073 /// # arguments 1074 /// 1075 /// * `vlength` - Vector length. Caller must ensure that number of 1076 /// quadrature points is a multiple of vlength. 1077 /// * `f` - Boxed closure to evaluate weak form at quadrature points. 1078 /// 1079 /// ``` 1080 /// # use libceed::prelude::*; 1081 /// # fn main() -> libceed::Result<()> { 1082 /// # let ceed = libceed::Ceed::default_init(); 1083 /// let mut user_f = |[u, weights, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1084 /// // Iterate over quadrature points 1085 /// v.iter_mut() 1086 /// .zip(u.iter().zip(weights.iter())) 1087 /// .for_each(|(v, (u, w))| *v = u * w); 1088 /// 1089 /// // Return clean error code 1090 /// 0 1091 /// }; 1092 /// 1093 /// let qf = ceed.q_function_interior(1, Box::new(user_f))?; 1094 /// # Ok(()) 1095 /// # } 1096 /// ``` 1097 pub fn q_function_interior<'a>( 1098 &self, 1099 vlength: usize, 1100 f: Box<qfunction::QFunctionUserClosure>, 1101 ) -> Result<QFunction<'a>> { 1102 QFunction::create(self, vlength, f) 1103 } 1104 1105 /// Returns a QFunction for evaluating interior (volumetric) terms 1106 /// created by name 1107 /// 1108 /// # arguments 1109 /// 1110 /// * `name` - name of QFunction from libCEED gallery 1111 /// 1112 /// ``` 1113 /// # use libceed::prelude::*; 1114 /// # fn main() -> libceed::Result<()> { 1115 /// # let ceed = libceed::Ceed::default_init(); 1116 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1117 /// # Ok(()) 1118 /// # } 1119 /// ``` 1120 pub fn q_function_interior_by_name<'a>(&self, name: &str) -> Result<QFunctionByName<'a>> { 1121 QFunctionByName::create(self, name) 1122 } 1123 1124 /// Returns an Operator and associate a QFunction. A Basis and 1125 /// ElemRestriction can be associated with QFunction fields via 1126 /// set_field(). 1127 /// 1128 /// # arguments 1129 /// 1130 /// * `qf` - QFunction defining the action of the operator at quadrature 1131 /// points 1132 /// * `dqf` - QFunction defining the action of the Jacobian of the qf (or 1133 /// qfunction_none) 1134 /// * `dqfT` - QFunction defining the action of the transpose of the 1135 /// Jacobian of the qf (or qfunction_none) 1136 /// 1137 /// ``` 1138 /// # use libceed::prelude::*; 1139 /// # fn main() -> libceed::Result<()> { 1140 /// # let ceed = libceed::Ceed::default_init(); 1141 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1142 /// let op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1143 /// # Ok(()) 1144 /// # } 1145 /// ``` 1146 pub fn operator<'a, 'b>( 1147 &self, 1148 qf: impl Into<QFunctionOpt<'b>>, 1149 dqf: impl Into<QFunctionOpt<'b>>, 1150 dqfT: impl Into<QFunctionOpt<'b>>, 1151 ) -> Result<Operator<'a>> { 1152 Operator::create(self, qf, dqf, dqfT) 1153 } 1154 1155 /// Returns an Operator that composes the action of several Operators 1156 /// 1157 /// ``` 1158 /// # use libceed::prelude::*; 1159 /// # fn main() -> libceed::Result<()> { 1160 /// # let ceed = libceed::Ceed::default_init(); 1161 /// let op = ceed.composite_operator()?; 1162 /// # Ok(()) 1163 /// # } 1164 /// ``` 1165 pub fn composite_operator<'a>(&self) -> Result<CompositeOperator<'a>> { 1166 CompositeOperator::create(self) 1167 } 1168 } 1169 1170 // ----------------------------------------------------------------------------- 1171 // Tests 1172 // ----------------------------------------------------------------------------- 1173 #[cfg(test)] 1174 mod tests { 1175 use super::*; 1176 1177 fn ceed_t501() -> Result<()> { 1178 let resource = "/cpu/self/ref/blocked"; 1179 let ceed = Ceed::init(resource); 1180 let nelem = 4; 1181 let p = 3; 1182 let q = 4; 1183 let ndofs = p * nelem - nelem + 1; 1184 1185 // Vectors 1186 let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1187 let mut qdata = ceed.vector(nelem * q)?; 1188 qdata.set_value(0.0)?; 1189 let mut u = ceed.vector(ndofs)?; 1190 u.set_value(1.0)?; 1191 let mut v = ceed.vector(ndofs)?; 1192 v.set_value(0.0)?; 1193 1194 // Restrictions 1195 let mut indx: Vec<i32> = vec![0; 2 * nelem]; 1196 for i in 0..nelem { 1197 for j in 0..2 { 1198 indx[2 * i + j] = (i + j) as i32; 1199 } 1200 } 1201 let rx = ceed.elem_restriction(nelem, 2, 1, 1, nelem + 1, MemType::Host, &indx)?; 1202 let mut indu: Vec<i32> = vec![0; p * nelem]; 1203 for i in 0..nelem { 1204 for j in 0..p { 1205 indu[p * i + j] = (i + j) as i32; 1206 } 1207 } 1208 let ru = ceed.elem_restriction(nelem, p, 1, 1, ndofs, MemType::Host, &indu)?; 1209 let strides: [i32; 3] = [1, q as i32, q as i32]; 1210 let rq = ceed.strided_elem_restriction(nelem, q, 1, q * nelem, strides)?; 1211 1212 // Bases 1213 let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1214 let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1215 1216 // Build quadrature data 1217 let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1218 ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1219 .field("dx", &rx, &bx, VectorOpt::Active)? 1220 .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1221 .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1222 .apply(&x, &mut qdata)?; 1223 1224 // Mass operator 1225 let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1226 let op_mass = ceed 1227 .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1228 .field("u", &ru, &bu, VectorOpt::Active)? 1229 .field("qdata", &rq, BasisOpt::None, &qdata)? 1230 .field("v", &ru, &bu, VectorOpt::Active)? 1231 .check()?; 1232 1233 v.set_value(0.0)?; 1234 op_mass.apply(&u, &mut v)?; 1235 1236 // Check 1237 let sum: Scalar = v.view()?.iter().sum(); 1238 let error: Scalar = (sum - 2.0).abs(); 1239 assert!( 1240 error < 50.0 * EPSILON, 1241 "Incorrect interval length computed. Expected: 2.0, Found: {}, Error: {:.12e}", 1242 sum, 1243 error 1244 ); 1245 Ok(()) 1246 } 1247 1248 #[test] 1249 fn test_ceed_t501() { 1250 assert!(ceed_t501().is_ok()); 1251 } 1252 } 1253 1254 // ----------------------------------------------------------------------------- 1255