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