1 // Copyright (c) 2017-2025, 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 //! A Ceed Operator defines the finite/spectral element operator associated to a 9 //! Ceed QFunction. A Ceed Operator connects Ceed ElemRestrictions, 10 //! Ceed Bases, and Ceed QFunctions. 11 12 use crate::{ 13 basis::{Basis, BasisOpt}, 14 elem_restriction::{ElemRestriction, ElemRestrictionOpt}, 15 prelude::*, 16 qfunction::QFunctionOpt, 17 vector::{Vector, VectorOpt}, 18 }; 19 20 // ----------------------------------------------------------------------------- 21 // Operator Field context wrapper 22 // ----------------------------------------------------------------------------- 23 #[derive(Debug)] 24 pub struct OperatorField<'a> { 25 pub(crate) ptr: bind_ceed::CeedOperatorField, 26 pub(crate) vector: crate::Vector<'a>, 27 pub(crate) elem_restriction: crate::ElemRestriction<'a>, 28 pub(crate) basis: crate::Basis<'a>, 29 _lifeline: PhantomData<&'a ()>, 30 } 31 32 // ----------------------------------------------------------------------------- 33 // Implementations 34 // ----------------------------------------------------------------------------- 35 impl<'a> OperatorField<'a> { 36 pub(crate) unsafe fn from_raw( 37 ptr: bind_ceed::CeedOperatorField, 38 ceed: crate::Ceed, 39 ) -> crate::Result<Self> { 40 let vector = { 41 let mut vector_ptr = std::ptr::null_mut(); 42 ceed.check_error(bind_ceed::CeedOperatorFieldGetVector(ptr, &mut vector_ptr))?; 43 crate::Vector::from_raw(vector_ptr)? 44 }; 45 let elem_restriction = { 46 let mut elem_restriction_ptr = std::ptr::null_mut(); 47 ceed.check_error(bind_ceed::CeedOperatorFieldGetElemRestriction( 48 ptr, 49 &mut elem_restriction_ptr, 50 ))?; 51 crate::ElemRestriction::from_raw(elem_restriction_ptr)? 52 }; 53 let basis = { 54 let mut basis_ptr = std::ptr::null_mut(); 55 ceed.check_error(bind_ceed::CeedOperatorFieldGetBasis(ptr, &mut basis_ptr))?; 56 crate::Basis::from_raw(basis_ptr)? 57 }; 58 Ok(Self { 59 ptr, 60 vector, 61 elem_restriction, 62 basis, 63 _lifeline: PhantomData, 64 }) 65 } 66 67 /// Get the name of an OperatorField 68 /// 69 /// ``` 70 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 71 /// # fn main() -> libceed::Result<()> { 72 /// # let ceed = libceed::Ceed::default_init(); 73 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 74 /// 75 /// // Operator field arguments 76 /// let ne = 3; 77 /// let q = 4_usize; 78 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 79 /// for i in 0..ne { 80 /// ind[2 * i + 0] = i as i32; 81 /// ind[2 * i + 1] = (i + 1) as i32; 82 /// } 83 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 84 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 85 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 86 /// 87 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 88 /// 89 /// // Operator fields 90 /// let op = ceed 91 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 92 /// .field("dx", &r, &b, VectorOpt::Active)? 93 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 94 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 95 /// 96 /// let inputs = op.inputs()?; 97 /// 98 /// assert_eq!(inputs[0].name(), "dx", "Incorrect input name"); 99 /// assert_eq!(inputs[1].name(), "weights", "Incorrect input name"); 100 /// # Ok(()) 101 /// # } 102 /// ``` 103 pub fn name(&self) -> &str { 104 let mut name_ptr: *mut std::os::raw::c_char = std::ptr::null_mut(); 105 unsafe { 106 bind_ceed::CeedOperatorFieldGetName( 107 self.ptr, 108 &mut name_ptr as *const _ as *mut *const _, 109 ); 110 } 111 unsafe { CStr::from_ptr(name_ptr) }.to_str().unwrap() 112 } 113 114 /// Get the ElemRestriction of an OperatorField 115 /// 116 /// ``` 117 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 118 /// # fn main() -> libceed::Result<()> { 119 /// # let ceed = libceed::Ceed::default_init(); 120 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 121 /// 122 /// // Operator field arguments 123 /// let ne = 3; 124 /// let q = 4_usize; 125 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 126 /// for i in 0..ne { 127 /// ind[2 * i + 0] = i as i32; 128 /// ind[2 * i + 1] = (i + 1) as i32; 129 /// } 130 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 131 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 132 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 133 /// 134 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 135 /// 136 /// // Operator fields 137 /// let op = ceed 138 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 139 /// .field("dx", &r, &b, VectorOpt::Active)? 140 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 141 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 142 /// 143 /// let inputs = op.inputs()?; 144 /// 145 /// assert!( 146 /// inputs[0].elem_restriction().is_some(), 147 /// "Incorrect field ElemRestriction" 148 /// ); 149 /// if let ElemRestrictionOpt::Some(r) = inputs[0].elem_restriction() { 150 /// assert_eq!( 151 /// r.num_elements(), 152 /// ne, 153 /// "Incorrect field ElemRestriction number of elements" 154 /// ); 155 /// } 156 /// 157 /// assert!( 158 /// inputs[1].elem_restriction().is_none(), 159 /// "Incorrect field ElemRestriction" 160 /// ); 161 /// 162 /// let outputs = op.outputs()?; 163 /// 164 /// assert!( 165 /// outputs[0].elem_restriction().is_some(), 166 /// "Incorrect field ElemRestriction" 167 /// ); 168 /// if let ElemRestrictionOpt::Some(r) = outputs[0].elem_restriction() { 169 /// assert_eq!( 170 /// r.num_elements(), 171 /// ne, 172 /// "Incorrect field ElemRestriction number of elements" 173 /// ); 174 /// } 175 /// # Ok(()) 176 /// # } 177 /// ``` 178 pub fn elem_restriction(&self) -> ElemRestrictionOpt { 179 if self.elem_restriction.ptr == unsafe { bind_ceed::CEED_ELEMRESTRICTION_NONE } { 180 ElemRestrictionOpt::None 181 } else { 182 ElemRestrictionOpt::Some(&self.elem_restriction) 183 } 184 } 185 186 /// Get the Basis of an OperatorField 187 /// 188 /// ``` 189 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 190 /// # fn main() -> libceed::Result<()> { 191 /// # let ceed = libceed::Ceed::default_init(); 192 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 193 /// 194 /// // Operator field arguments 195 /// let ne = 3; 196 /// let q = 4_usize; 197 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 198 /// for i in 0..ne { 199 /// ind[2 * i + 0] = i as i32; 200 /// ind[2 * i + 1] = (i + 1) as i32; 201 /// } 202 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 203 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 204 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 205 /// 206 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 207 /// 208 /// // Operator fields 209 /// let op = ceed 210 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 211 /// .field("dx", &r, &b, VectorOpt::Active)? 212 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 213 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 214 /// 215 /// let inputs = op.inputs()?; 216 /// 217 /// assert!(inputs[0].basis().is_some(), "Incorrect field Basis"); 218 /// if let BasisOpt::Some(b) = inputs[0].basis() { 219 /// assert_eq!( 220 /// b.num_quadrature_points(), 221 /// q, 222 /// "Incorrect field Basis number of quadrature points" 223 /// ); 224 /// } 225 /// assert!(inputs[1].basis().is_some(), "Incorrect field Basis"); 226 /// if let BasisOpt::Some(b) = inputs[1].basis() { 227 /// assert_eq!( 228 /// b.num_quadrature_points(), 229 /// q, 230 /// "Incorrect field Basis number of quadrature points" 231 /// ); 232 /// } 233 /// 234 /// let outputs = op.outputs()?; 235 /// 236 /// assert!(outputs[0].basis().is_none(), "Incorrect field Basis"); 237 /// # Ok(()) 238 /// # } 239 /// ``` 240 pub fn basis(&self) -> BasisOpt { 241 if self.basis.ptr == unsafe { bind_ceed::CEED_BASIS_NONE } { 242 BasisOpt::None 243 } else { 244 BasisOpt::Some(&self.basis) 245 } 246 } 247 248 /// Get the Vector of an OperatorField 249 /// 250 /// ``` 251 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 252 /// # fn main() -> libceed::Result<()> { 253 /// # let ceed = libceed::Ceed::default_init(); 254 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 255 /// 256 /// // Operator field arguments 257 /// let ne = 3; 258 /// let q = 4_usize; 259 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 260 /// for i in 0..ne { 261 /// ind[2 * i + 0] = i as i32; 262 /// ind[2 * i + 1] = (i + 1) as i32; 263 /// } 264 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 265 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 266 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 267 /// 268 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 269 /// 270 /// // Operator fields 271 /// let op = ceed 272 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 273 /// .field("dx", &r, &b, VectorOpt::Active)? 274 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 275 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 276 /// 277 /// let inputs = op.inputs()?; 278 /// 279 /// assert!(inputs[0].vector().is_active(), "Incorrect field Vector"); 280 /// assert!(inputs[1].vector().is_none(), "Incorrect field Vector"); 281 /// 282 /// let outputs = op.outputs()?; 283 /// 284 /// assert!(outputs[0].vector().is_active(), "Incorrect field Vector"); 285 /// # Ok(()) 286 /// # } 287 /// ``` 288 pub fn vector(&self) -> VectorOpt { 289 if self.vector.ptr == unsafe { bind_ceed::CEED_VECTOR_ACTIVE } { 290 VectorOpt::Active 291 } else if self.vector.ptr == unsafe { bind_ceed::CEED_VECTOR_NONE } { 292 VectorOpt::None 293 } else { 294 VectorOpt::Some(&self.vector) 295 } 296 } 297 } 298 299 // ----------------------------------------------------------------------------- 300 // Operator context wrapper 301 // ----------------------------------------------------------------------------- 302 #[derive(Debug)] 303 pub(crate) struct OperatorCore<'a> { 304 ptr: bind_ceed::CeedOperator, 305 _lifeline: PhantomData<&'a ()>, 306 } 307 308 #[derive(Debug)] 309 pub struct Operator<'a> { 310 op_core: OperatorCore<'a>, 311 } 312 313 #[derive(Debug)] 314 pub struct CompositeOperator<'a> { 315 op_core: OperatorCore<'a>, 316 } 317 318 // ----------------------------------------------------------------------------- 319 // Destructor 320 // ----------------------------------------------------------------------------- 321 impl<'a> Drop for OperatorCore<'a> { 322 fn drop(&mut self) { 323 unsafe { 324 bind_ceed::CeedOperatorDestroy(&mut self.ptr); 325 } 326 } 327 } 328 329 // ----------------------------------------------------------------------------- 330 // Display 331 // ----------------------------------------------------------------------------- 332 impl<'a> fmt::Display for OperatorCore<'a> { 333 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 334 let mut ptr = std::ptr::null_mut(); 335 let mut sizeloc = crate::MAX_BUFFER_LENGTH; 336 let cstring = unsafe { 337 let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc); 338 bind_ceed::CeedOperatorView(self.ptr, file); 339 bind_ceed::fclose(file); 340 CString::from_raw(ptr) 341 }; 342 cstring.to_string_lossy().fmt(f) 343 } 344 } 345 346 /// View an Operator 347 /// 348 /// ``` 349 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 350 /// # fn main() -> libceed::Result<()> { 351 /// # let ceed = libceed::Ceed::default_init(); 352 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 353 /// 354 /// // Operator field arguments 355 /// let ne = 3; 356 /// let q = 4_usize; 357 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 358 /// for i in 0..ne { 359 /// ind[2 * i + 0] = i as i32; 360 /// ind[2 * i + 1] = (i + 1) as i32; 361 /// } 362 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 363 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 364 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 365 /// 366 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 367 /// 368 /// // Operator fields 369 /// let op = ceed 370 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 371 /// .name("mass")? 372 /// .field("dx", &r, &b, VectorOpt::Active)? 373 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 374 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 375 /// 376 /// println!("{}", op); 377 /// # Ok(()) 378 /// # } 379 /// ``` 380 impl<'a> fmt::Display for Operator<'a> { 381 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 382 self.op_core.fmt(f) 383 } 384 } 385 386 /// View a composite Operator 387 /// 388 /// ``` 389 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 390 /// # fn main() -> libceed::Result<()> { 391 /// # let ceed = libceed::Ceed::default_init(); 392 /// 393 /// // Sub operator field arguments 394 /// let ne = 3; 395 /// let q = 4_usize; 396 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 397 /// for i in 0..ne { 398 /// ind[2 * i + 0] = i as i32; 399 /// ind[2 * i + 1] = (i + 1) as i32; 400 /// } 401 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 402 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 403 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 404 /// 405 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 406 /// 407 /// let qdata_mass = ceed.vector(q * ne)?; 408 /// let qdata_diff = ceed.vector(q * ne)?; 409 /// 410 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 411 /// let op_mass = ceed 412 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 413 /// .name("Mass term")? 414 /// .field("u", &r, &b, VectorOpt::Active)? 415 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 416 /// .field("v", &r, &b, VectorOpt::Active)?; 417 /// 418 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 419 /// let op_diff = ceed 420 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 421 /// .name("Poisson term")? 422 /// .field("du", &r, &b, VectorOpt::Active)? 423 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 424 /// .field("dv", &r, &b, VectorOpt::Active)?; 425 /// 426 /// let op = ceed 427 /// .composite_operator()? 428 /// .name("Screened Poisson")? 429 /// .sub_operator(&op_mass)? 430 /// .sub_operator(&op_diff)?; 431 /// 432 /// println!("{}", op); 433 /// # Ok(()) 434 /// # } 435 /// ``` 436 impl<'a> fmt::Display for CompositeOperator<'a> { 437 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { 438 self.op_core.fmt(f) 439 } 440 } 441 442 // ----------------------------------------------------------------------------- 443 // Core functionality 444 // ----------------------------------------------------------------------------- 445 impl<'a> OperatorCore<'a> { 446 // Raw Ceed for error handling 447 #[doc(hidden)] 448 fn ceed(&self) -> bind_ceed::Ceed { 449 unsafe { bind_ceed::CeedOperatorReturnCeed(self.ptr) } 450 } 451 452 // Error handling 453 #[doc(hidden)] 454 fn check_error(&self, ierr: i32) -> crate::Result<i32> { 455 crate::check_error(|| self.ceed(), ierr) 456 } 457 458 // Common implementations 459 pub fn check(&self) -> crate::Result<i32> { 460 self.check_error(unsafe { bind_ceed::CeedOperatorCheckReady(self.ptr) }) 461 } 462 463 pub fn name(&self, name: &str) -> crate::Result<i32> { 464 let name_c = CString::new(name).expect("CString::new failed"); 465 self.check_error(unsafe { bind_ceed::CeedOperatorSetName(self.ptr, name_c.as_ptr()) }) 466 } 467 468 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 469 self.check_error(unsafe { 470 bind_ceed::CeedOperatorApply( 471 self.ptr, 472 input.ptr, 473 output.ptr, 474 bind_ceed::CEED_REQUEST_IMMEDIATE, 475 ) 476 }) 477 } 478 479 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 480 self.check_error(unsafe { 481 bind_ceed::CeedOperatorApplyAdd( 482 self.ptr, 483 input.ptr, 484 output.ptr, 485 bind_ceed::CEED_REQUEST_IMMEDIATE, 486 ) 487 }) 488 } 489 490 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 491 self.check_error(unsafe { 492 bind_ceed::CeedOperatorLinearAssembleDiagonal( 493 self.ptr, 494 assembled.ptr, 495 bind_ceed::CEED_REQUEST_IMMEDIATE, 496 ) 497 }) 498 } 499 500 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 501 self.check_error(unsafe { 502 bind_ceed::CeedOperatorLinearAssembleAddDiagonal( 503 self.ptr, 504 assembled.ptr, 505 bind_ceed::CEED_REQUEST_IMMEDIATE, 506 ) 507 }) 508 } 509 510 pub fn linear_assemble_point_block_diagonal( 511 &self, 512 assembled: &mut Vector, 513 ) -> crate::Result<i32> { 514 self.check_error(unsafe { 515 bind_ceed::CeedOperatorLinearAssemblePointBlockDiagonal( 516 self.ptr, 517 assembled.ptr, 518 bind_ceed::CEED_REQUEST_IMMEDIATE, 519 ) 520 }) 521 } 522 523 pub fn linear_assemble_add_point_block_diagonal( 524 &self, 525 assembled: &mut Vector, 526 ) -> crate::Result<i32> { 527 self.check_error(unsafe { 528 bind_ceed::CeedOperatorLinearAssembleAddPointBlockDiagonal( 529 self.ptr, 530 assembled.ptr, 531 bind_ceed::CEED_REQUEST_IMMEDIATE, 532 ) 533 }) 534 } 535 } 536 537 // ----------------------------------------------------------------------------- 538 // Operator 539 // ----------------------------------------------------------------------------- 540 impl<'a> Operator<'a> { 541 // Constructor 542 pub fn create<'b>( 543 ceed: &crate::Ceed, 544 qf: impl Into<QFunctionOpt<'b>>, 545 dqf: impl Into<QFunctionOpt<'b>>, 546 dqfT: impl Into<QFunctionOpt<'b>>, 547 ) -> crate::Result<Self> { 548 let mut ptr = std::ptr::null_mut(); 549 ceed.check_error(unsafe { 550 bind_ceed::CeedOperatorCreate( 551 ceed.ptr, 552 qf.into().to_raw(), 553 dqf.into().to_raw(), 554 dqfT.into().to_raw(), 555 &mut ptr, 556 ) 557 })?; 558 Ok(Self { 559 op_core: OperatorCore { 560 ptr, 561 _lifeline: PhantomData, 562 }, 563 }) 564 } 565 566 unsafe fn from_raw(ptr: bind_ceed::CeedOperator) -> crate::Result<Self> { 567 Ok(Self { 568 op_core: OperatorCore { 569 ptr, 570 _lifeline: PhantomData, 571 }, 572 }) 573 } 574 575 /// Set name for Operator printing 576 /// 577 /// * 'name' - Name to set 578 /// 579 /// ``` 580 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 581 /// # fn main() -> libceed::Result<()> { 582 /// # let ceed = libceed::Ceed::default_init(); 583 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 584 /// 585 /// // Operator field arguments 586 /// let ne = 3; 587 /// let q = 4_usize; 588 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 589 /// for i in 0..ne { 590 /// ind[2 * i + 0] = i as i32; 591 /// ind[2 * i + 1] = (i + 1) as i32; 592 /// } 593 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 594 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 595 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 596 /// 597 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 598 /// 599 /// // Operator fields 600 /// let op = ceed 601 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 602 /// .name("mass")? 603 /// .field("dx", &r, &b, VectorOpt::Active)? 604 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 605 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 606 /// # Ok(()) 607 /// # } 608 /// ``` 609 #[allow(unused_mut)] 610 pub fn name(mut self, name: &str) -> crate::Result<Self> { 611 self.op_core.name(name)?; 612 Ok(self) 613 } 614 615 /// Apply Operator to a vector 616 /// 617 /// * `input` - Input Vector 618 /// * `output` - Output Vector 619 /// 620 /// ``` 621 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 622 /// # fn main() -> libceed::Result<()> { 623 /// # let ceed = libceed::Ceed::default_init(); 624 /// let ne = 4; 625 /// let p = 3; 626 /// let q = 4; 627 /// let ndofs = p * ne - ne + 1; 628 /// 629 /// // Vectors 630 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 631 /// let mut qdata = ceed.vector(ne * q)?; 632 /// qdata.set_value(0.0); 633 /// let u = ceed.vector_from_slice(&vec![1.0; ndofs])?; 634 /// let mut v = ceed.vector(ndofs)?; 635 /// v.set_value(0.0); 636 /// 637 /// // Restrictions 638 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 639 /// for i in 0..ne { 640 /// indx[2 * i + 0] = i as i32; 641 /// indx[2 * i + 1] = (i + 1) as i32; 642 /// } 643 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 644 /// let mut indu: Vec<i32> = vec![0; p * ne]; 645 /// for i in 0..ne { 646 /// indu[p * i + 0] = i as i32; 647 /// indu[p * i + 1] = (i + 1) as i32; 648 /// indu[p * i + 2] = (i + 2) as i32; 649 /// } 650 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 651 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 652 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 653 /// 654 /// // Bases 655 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 656 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 657 /// 658 /// // Build quadrature data 659 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 660 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 661 /// .field("dx", &rx, &bx, VectorOpt::Active)? 662 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 663 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 664 /// .apply(&x, &mut qdata)?; 665 /// 666 /// // Mass operator 667 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 668 /// let op_mass = ceed 669 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 670 /// .field("u", &ru, &bu, VectorOpt::Active)? 671 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 672 /// .field("v", &ru, &bu, VectorOpt::Active)?; 673 /// 674 /// v.set_value(0.0); 675 /// op_mass.apply(&u, &mut v)?; 676 /// 677 /// // Check 678 /// let sum: Scalar = v.view()?.iter().sum(); 679 /// let error: Scalar = (sum - 2.0).abs(); 680 /// assert!( 681 /// error < 50.0 * libceed::EPSILON, 682 /// "Incorrect interval length computed. Expected: 2.0, Found: {}, Error: {:.12e}", 683 /// sum, 684 /// error 685 /// ); 686 /// # Ok(()) 687 /// # } 688 /// ``` 689 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 690 self.op_core.apply(input, output) 691 } 692 693 /// Apply Operator to a vector and add result to output vector 694 /// 695 /// * `input` - Input Vector 696 /// * `output` - Output Vector 697 /// 698 /// ``` 699 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 700 /// # fn main() -> libceed::Result<()> { 701 /// # let ceed = libceed::Ceed::default_init(); 702 /// let ne = 4; 703 /// let p = 3; 704 /// let q = 4; 705 /// let ndofs = p * ne - ne + 1; 706 /// 707 /// // Vectors 708 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 709 /// let mut qdata = ceed.vector(ne * q)?; 710 /// qdata.set_value(0.0); 711 /// let u = ceed.vector_from_slice(&vec![1.0; ndofs])?; 712 /// let mut v = ceed.vector(ndofs)?; 713 /// 714 /// // Restrictions 715 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 716 /// for i in 0..ne { 717 /// indx[2 * i + 0] = i as i32; 718 /// indx[2 * i + 1] = (i + 1) as i32; 719 /// } 720 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 721 /// let mut indu: Vec<i32> = vec![0; p * ne]; 722 /// for i in 0..ne { 723 /// indu[p * i + 0] = i as i32; 724 /// indu[p * i + 1] = (i + 1) as i32; 725 /// indu[p * i + 2] = (i + 2) as i32; 726 /// } 727 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 728 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 729 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 730 /// 731 /// // Bases 732 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 733 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 734 /// 735 /// // Build quadrature data 736 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 737 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 738 /// .field("dx", &rx, &bx, VectorOpt::Active)? 739 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 740 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 741 /// .apply(&x, &mut qdata)?; 742 /// 743 /// // Mass operator 744 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 745 /// let op_mass = ceed 746 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 747 /// .field("u", &ru, &bu, VectorOpt::Active)? 748 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 749 /// .field("v", &ru, &bu, VectorOpt::Active)?; 750 /// 751 /// v.set_value(1.0); 752 /// op_mass.apply_add(&u, &mut v)?; 753 /// 754 /// // Check 755 /// let sum: Scalar = v.view()?.iter().sum(); 756 /// assert!( 757 /// (sum - (2.0 + ndofs as Scalar)).abs() < 10.0 * libceed::EPSILON, 758 /// "Incorrect interval length computed and added" 759 /// ); 760 /// # Ok(()) 761 /// # } 762 /// ``` 763 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 764 self.op_core.apply_add(input, output) 765 } 766 767 /// Provide a field to a Operator for use by its QFunction 768 /// 769 /// * `fieldname` - Name of the field (to be matched with the name used by 770 /// the QFunction) 771 /// * `r` - ElemRestriction 772 /// * `b` - Basis in which the field resides or `BasisOpt::None` 773 /// * `v` - Vector to be used by Operator or `VectorOpt::Active` if field 774 /// is active or `VectorOpt::None` if using `Weight` with the 775 /// QFunction 776 /// 777 /// 778 /// ``` 779 /// # use libceed::{prelude::*, MemType, QFunctionOpt, QuadMode, VectorOpt}; 780 /// # fn main() -> libceed::Result<()> { 781 /// # let ceed = libceed::Ceed::default_init(); 782 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 783 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 784 /// 785 /// // Operator field arguments 786 /// let ne = 3; 787 /// let q = 4; 788 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 789 /// for i in 0..ne { 790 /// ind[2 * i + 0] = i as i32; 791 /// ind[2 * i + 1] = (i + 1) as i32; 792 /// } 793 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 794 /// 795 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 796 /// 797 /// // Operator field 798 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 799 /// # Ok(()) 800 /// # } 801 /// ``` 802 #[allow(unused_mut)] 803 pub fn field<'b>( 804 mut self, 805 fieldname: &str, 806 r: impl Into<ElemRestrictionOpt<'b>>, 807 b: impl Into<BasisOpt<'b>>, 808 v: impl Into<VectorOpt<'b>>, 809 ) -> crate::Result<Self> { 810 let fieldname = CString::new(fieldname).expect("CString::new failed"); 811 let fieldname = fieldname.as_ptr(); 812 self.op_core.check_error(unsafe { 813 bind_ceed::CeedOperatorSetField( 814 self.op_core.ptr, 815 fieldname, 816 r.into().to_raw(), 817 b.into().to_raw(), 818 v.into().to_raw(), 819 ) 820 })?; 821 Ok(self) 822 } 823 824 /// Get a slice of Operator inputs 825 /// 826 /// ``` 827 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 828 /// # fn main() -> libceed::Result<()> { 829 /// # let ceed = libceed::Ceed::default_init(); 830 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 831 /// 832 /// // Operator field arguments 833 /// let ne = 3; 834 /// let q = 4_usize; 835 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 836 /// for i in 0..ne { 837 /// ind[2 * i + 0] = i as i32; 838 /// ind[2 * i + 1] = (i + 1) as i32; 839 /// } 840 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 841 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 842 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 843 /// 844 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 845 /// 846 /// // Operator fields 847 /// let op = ceed 848 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 849 /// .field("dx", &r, &b, VectorOpt::Active)? 850 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 851 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 852 /// 853 /// let inputs = op.inputs()?; 854 /// 855 /// assert_eq!(inputs.len(), 2, "Incorrect inputs array"); 856 /// # Ok(()) 857 /// # } 858 /// ``` 859 pub fn inputs(&self) -> crate::Result<Vec<crate::OperatorField>> { 860 // Get array of raw C pointers for inputs 861 let mut num_inputs = 0; 862 let mut inputs_ptr = std::ptr::null_mut(); 863 self.op_core.check_error(unsafe { 864 bind_ceed::CeedOperatorGetFields( 865 self.op_core.ptr, 866 &mut num_inputs, 867 &mut inputs_ptr, 868 std::ptr::null_mut() as *mut bind_ceed::CeedInt, 869 std::ptr::null_mut(), 870 ) 871 })?; 872 // Convert raw C pointers to fixed length slice 873 let inputs_slice = unsafe { 874 std::slice::from_raw_parts( 875 inputs_ptr as *mut bind_ceed::CeedOperatorField, 876 num_inputs as usize, 877 ) 878 }; 879 // And finally build vec 880 let ceed = { 881 let ceed_raw = self.op_core.ceed(); 882 let mut ptr = std::ptr::null_mut(); 883 unsafe { 884 bind_ceed::CeedReferenceCopy(ceed_raw, &mut ptr); // refcount 885 } 886 crate::Ceed { ptr } 887 }; 888 let inputs = (0..num_inputs as usize) 889 .map(|i| unsafe { crate::OperatorField::from_raw(inputs_slice[i], ceed.clone()) }) 890 .collect::<crate::Result<Vec<_>>>()?; 891 Ok(inputs) 892 } 893 894 /// Get a slice of Operator outputs 895 /// 896 /// ``` 897 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 898 /// # fn main() -> libceed::Result<()> { 899 /// # let ceed = libceed::Ceed::default_init(); 900 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 901 /// 902 /// // Operator field arguments 903 /// let ne = 3; 904 /// let q = 4_usize; 905 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 906 /// for i in 0..ne { 907 /// ind[2 * i + 0] = i as i32; 908 /// ind[2 * i + 1] = (i + 1) as i32; 909 /// } 910 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 911 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 912 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 913 /// 914 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 915 /// 916 /// // Operator fields 917 /// let op = ceed 918 /// .operator(&qf, QFunctionOpt::None, QFunctionOpt::None)? 919 /// .field("dx", &r, &b, VectorOpt::Active)? 920 /// .field("weights", ElemRestrictionOpt::None, &b, VectorOpt::None)? 921 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 922 /// 923 /// let outputs = op.outputs()?; 924 /// 925 /// assert_eq!(outputs.len(), 1, "Incorrect outputs array"); 926 /// # Ok(()) 927 /// # } 928 /// ``` 929 pub fn outputs(&self) -> crate::Result<Vec<crate::OperatorField>> { 930 // Get array of raw C pointers for outputs 931 let mut num_outputs = 0; 932 let mut outputs_ptr = std::ptr::null_mut(); 933 self.op_core.check_error(unsafe { 934 bind_ceed::CeedOperatorGetFields( 935 self.op_core.ptr, 936 std::ptr::null_mut() as *mut bind_ceed::CeedInt, 937 std::ptr::null_mut(), 938 &mut num_outputs, 939 &mut outputs_ptr, 940 ) 941 })?; 942 // Convert raw C pointers to fixed length slice 943 let outputs_slice = unsafe { 944 std::slice::from_raw_parts( 945 outputs_ptr as *mut bind_ceed::CeedOperatorField, 946 num_outputs as usize, 947 ) 948 }; 949 // And finally build vec 950 let ceed = { 951 let ceed_raw = self.op_core.ceed(); 952 let mut ptr = std::ptr::null_mut(); 953 unsafe { 954 bind_ceed::CeedReferenceCopy(ceed_raw, &mut ptr); // refcount 955 } 956 crate::Ceed { ptr } 957 }; 958 let outputs = (0..num_outputs as usize) 959 .map(|i| unsafe { crate::OperatorField::from_raw(outputs_slice[i], ceed.clone()) }) 960 .collect::<crate::Result<Vec<_>>>()?; 961 Ok(outputs) 962 } 963 964 /// Check if Operator is setup correctly 965 /// 966 /// ``` 967 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 968 /// # fn main() -> libceed::Result<()> { 969 /// # let ceed = libceed::Ceed::default_init(); 970 /// let ne = 4; 971 /// let p = 3; 972 /// let q = 4; 973 /// let ndofs = p * ne - ne + 1; 974 /// 975 /// // Restrictions 976 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 977 /// for i in 0..ne { 978 /// indx[2 * i + 0] = i as i32; 979 /// indx[2 * i + 1] = (i + 1) as i32; 980 /// } 981 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 982 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 983 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 984 /// 985 /// // Bases 986 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 987 /// 988 /// // Build quadrature data 989 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 990 /// let op_build = ceed 991 /// .operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 992 /// .field("dx", &rx, &bx, VectorOpt::Active)? 993 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 994 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 995 /// 996 /// // Check 997 /// assert!(op_build.check().is_ok(), "Operator setup incorrectly"); 998 /// # Ok(()) 999 /// # } 1000 /// ``` 1001 pub fn check(self) -> crate::Result<Self> { 1002 self.op_core.check()?; 1003 Ok(self) 1004 } 1005 1006 /// Get the number of elements associated with an Operator 1007 /// 1008 /// 1009 /// ``` 1010 /// # use libceed::{prelude::*, MemType, QFunctionOpt, QuadMode, VectorOpt}; 1011 /// # fn main() -> libceed::Result<()> { 1012 /// # let ceed = libceed::Ceed::default_init(); 1013 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1014 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1015 /// 1016 /// // Operator field arguments 1017 /// let ne = 3; 1018 /// let q = 4; 1019 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 1020 /// for i in 0..ne { 1021 /// ind[2 * i + 0] = i as i32; 1022 /// ind[2 * i + 1] = (i + 1) as i32; 1023 /// } 1024 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 1025 /// 1026 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1027 /// 1028 /// // Operator field 1029 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 1030 /// 1031 /// // Check number of elements 1032 /// assert_eq!(op.num_elements(), ne, "Incorrect number of elements"); 1033 /// # Ok(()) 1034 /// # } 1035 /// ``` 1036 pub fn num_elements(&self) -> usize { 1037 let mut nelem = 0; 1038 unsafe { bind_ceed::CeedOperatorGetNumElements(self.op_core.ptr, &mut nelem) }; 1039 usize::try_from(nelem).unwrap() 1040 } 1041 1042 /// Get the number of quadrature points associated with each element of 1043 /// an Operator 1044 /// 1045 /// 1046 /// ``` 1047 /// # use libceed::{prelude::*, MemType, QFunctionOpt, QuadMode, VectorOpt}; 1048 /// # fn main() -> libceed::Result<()> { 1049 /// # let ceed = libceed::Ceed::default_init(); 1050 /// let qf = ceed.q_function_interior_by_name("Mass1DBuild")?; 1051 /// let mut op = ceed.operator(&qf, QFunctionOpt::None, QFunctionOpt::None)?; 1052 /// 1053 /// // Operator field arguments 1054 /// let ne = 3; 1055 /// let q = 4; 1056 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 1057 /// for i in 0..ne { 1058 /// ind[2 * i + 0] = i as i32; 1059 /// ind[2 * i + 1] = (i + 1) as i32; 1060 /// } 1061 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 1062 /// 1063 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1064 /// 1065 /// // Operator field 1066 /// op = op.field("dx", &r, &b, VectorOpt::Active)?; 1067 /// 1068 /// // Check number of quadrature points 1069 /// assert_eq!( 1070 /// op.num_quadrature_points(), 1071 /// q, 1072 /// "Incorrect number of quadrature points" 1073 /// ); 1074 /// # Ok(()) 1075 /// # } 1076 /// ``` 1077 pub fn num_quadrature_points(&self) -> usize { 1078 let mut Q = 0; 1079 unsafe { bind_ceed::CeedOperatorGetNumQuadraturePoints(self.op_core.ptr, &mut Q) }; 1080 usize::try_from(Q).unwrap() 1081 } 1082 1083 /// Assemble the diagonal of a square linear Operator 1084 /// 1085 /// This overwrites a Vector with the diagonal of a linear Operator. 1086 /// 1087 /// Note: Currently only non-composite Operators with a single field and 1088 /// composite Operators with single field sub-operators are supported. 1089 /// 1090 /// * `op` - Operator to assemble QFunction 1091 /// * `assembled` - Vector to store assembled Operator diagonal 1092 /// 1093 /// ``` 1094 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 1095 /// # fn main() -> libceed::Result<()> { 1096 /// # let ceed = libceed::Ceed::default_init(); 1097 /// let ne = 4; 1098 /// let p = 3; 1099 /// let q = 4; 1100 /// let ndofs = p * ne - ne + 1; 1101 /// 1102 /// // Vectors 1103 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1104 /// let mut qdata = ceed.vector(ne * q)?; 1105 /// qdata.set_value(0.0); 1106 /// let mut u = ceed.vector(ndofs)?; 1107 /// u.set_value(1.0); 1108 /// let mut v = ceed.vector(ndofs)?; 1109 /// v.set_value(0.0); 1110 /// 1111 /// // Restrictions 1112 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1113 /// for i in 0..ne { 1114 /// indx[2 * i + 0] = i as i32; 1115 /// indx[2 * i + 1] = (i + 1) as i32; 1116 /// } 1117 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1118 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1119 /// for i in 0..ne { 1120 /// indu[p * i + 0] = (2 * i) as i32; 1121 /// indu[p * i + 1] = (2 * i + 1) as i32; 1122 /// indu[p * i + 2] = (2 * i + 2) as i32; 1123 /// } 1124 /// let ru = ceed.elem_restriction(ne, p, 1, 1, ndofs, MemType::Host, &indu)?; 1125 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1126 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1127 /// 1128 /// // Bases 1129 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1130 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1131 /// 1132 /// // Build quadrature data 1133 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1134 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1135 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1136 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1137 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1138 /// .apply(&x, &mut qdata)?; 1139 /// 1140 /// // Mass operator 1141 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1142 /// let op_mass = ceed 1143 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1144 /// .field("u", &ru, &bu, VectorOpt::Active)? 1145 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1146 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1147 /// 1148 /// // Diagonal 1149 /// let mut diag = ceed.vector(ndofs)?; 1150 /// diag.set_value(0.0); 1151 /// op_mass.linear_assemble_diagonal(&mut diag)?; 1152 /// 1153 /// // Manual diagonal computation 1154 /// let mut true_diag = ceed.vector(ndofs)?; 1155 /// true_diag.set_value(0.0)?; 1156 /// for i in 0..ndofs { 1157 /// u.set_value(0.0); 1158 /// { 1159 /// let mut u_array = u.view_mut()?; 1160 /// u_array[i] = 1.; 1161 /// } 1162 /// 1163 /// op_mass.apply(&u, &mut v)?; 1164 /// 1165 /// { 1166 /// let v_array = v.view()?; 1167 /// let mut true_array = true_diag.view_mut()?; 1168 /// true_array[i] = v_array[i]; 1169 /// } 1170 /// } 1171 /// 1172 /// // Check 1173 /// diag.view()? 1174 /// .iter() 1175 /// .zip(true_diag.view()?.iter()) 1176 /// .for_each(|(computed, actual)| { 1177 /// assert!( 1178 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1179 /// "Diagonal entry incorrect" 1180 /// ); 1181 /// }); 1182 /// # Ok(()) 1183 /// # } 1184 /// ``` 1185 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 1186 self.op_core.linear_assemble_diagonal(assembled) 1187 } 1188 1189 /// Assemble the diagonal of a square linear Operator 1190 /// 1191 /// This sums into a Vector with the diagonal of a linear Operator. 1192 /// 1193 /// Note: Currently only non-composite Operators with a single field and 1194 /// composite Operators with single field sub-operators are supported. 1195 /// 1196 /// * `op` - Operator to assemble QFunction 1197 /// * `assembled` - Vector to store assembled Operator diagonal 1198 /// 1199 /// 1200 /// ``` 1201 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 1202 /// # fn main() -> libceed::Result<()> { 1203 /// # let ceed = libceed::Ceed::default_init(); 1204 /// let ne = 4; 1205 /// let p = 3; 1206 /// let q = 4; 1207 /// let ndofs = p * ne - ne + 1; 1208 /// 1209 /// // Vectors 1210 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1211 /// let mut qdata = ceed.vector(ne * q)?; 1212 /// qdata.set_value(0.0); 1213 /// let mut u = ceed.vector(ndofs)?; 1214 /// u.set_value(1.0); 1215 /// let mut v = ceed.vector(ndofs)?; 1216 /// v.set_value(0.0); 1217 /// 1218 /// // Restrictions 1219 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1220 /// for i in 0..ne { 1221 /// indx[2 * i + 0] = i as i32; 1222 /// indx[2 * i + 1] = (i + 1) as i32; 1223 /// } 1224 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1225 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1226 /// for i in 0..ne { 1227 /// indu[p * i + 0] = (2 * i) as i32; 1228 /// indu[p * i + 1] = (2 * i + 1) as i32; 1229 /// indu[p * i + 2] = (2 * i + 2) as i32; 1230 /// } 1231 /// let ru = ceed.elem_restriction(ne, p, 1, 1, ndofs, MemType::Host, &indu)?; 1232 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1233 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1234 /// 1235 /// // Bases 1236 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1237 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 1238 /// 1239 /// // Build quadrature data 1240 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1241 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1242 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1243 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1244 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1245 /// .apply(&x, &mut qdata)?; 1246 /// 1247 /// // Mass operator 1248 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1249 /// let op_mass = ceed 1250 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1251 /// .field("u", &ru, &bu, VectorOpt::Active)? 1252 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1253 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1254 /// 1255 /// // Diagonal 1256 /// let mut diag = ceed.vector(ndofs)?; 1257 /// diag.set_value(1.0); 1258 /// op_mass.linear_assemble_add_diagonal(&mut diag)?; 1259 /// 1260 /// // Manual diagonal computation 1261 /// let mut true_diag = ceed.vector(ndofs)?; 1262 /// true_diag.set_value(0.0)?; 1263 /// for i in 0..ndofs { 1264 /// u.set_value(0.0); 1265 /// { 1266 /// let mut u_array = u.view_mut()?; 1267 /// u_array[i] = 1.; 1268 /// } 1269 /// 1270 /// op_mass.apply(&u, &mut v)?; 1271 /// 1272 /// { 1273 /// let v_array = v.view()?; 1274 /// let mut true_array = true_diag.view_mut()?; 1275 /// true_array[i] = v_array[i] + 1.0; 1276 /// } 1277 /// } 1278 /// 1279 /// // Check 1280 /// diag.view()? 1281 /// .iter() 1282 /// .zip(true_diag.view()?.iter()) 1283 /// .for_each(|(computed, actual)| { 1284 /// assert!( 1285 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1286 /// "Diagonal entry incorrect" 1287 /// ); 1288 /// }); 1289 /// # Ok(()) 1290 /// # } 1291 /// ``` 1292 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 1293 self.op_core.linear_assemble_add_diagonal(assembled) 1294 } 1295 1296 /// Assemble the point block diagonal of a square linear Operator 1297 /// 1298 /// This overwrites a Vector with the point block diagonal of a linear 1299 /// Operator. 1300 /// 1301 /// Note: Currently only non-composite Operators with a single field and 1302 /// composite Operators with single field sub-operators are supported. 1303 /// 1304 /// * `op` - Operator to assemble QFunction 1305 /// * `assembled` - Vector to store assembled CeedOperator point block 1306 /// diagonal, provided in row-major form with an 1307 /// `ncomp * ncomp` block at each node. The dimensions of 1308 /// this vector are derived from the active vector for 1309 /// the CeedOperator. The array has shape 1310 /// `[nodes, component out, component in]`. 1311 /// 1312 /// ``` 1313 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, EvalMode, MemType, QFunctionInputs, QFunctionOpt, QFunctionOutputs, QuadMode, VectorOpt}; 1314 /// # fn main() -> libceed::Result<()> { 1315 /// # let ceed = libceed::Ceed::default_init(); 1316 /// let ne = 4; 1317 /// let p = 3; 1318 /// let q = 4; 1319 /// let ncomp = 2; 1320 /// let ndofs = p * ne - ne + 1; 1321 /// 1322 /// // Vectors 1323 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1324 /// let mut qdata = ceed.vector(ne * q)?; 1325 /// qdata.set_value(0.0); 1326 /// let mut u = ceed.vector(ncomp * ndofs)?; 1327 /// u.set_value(1.0); 1328 /// let mut v = ceed.vector(ncomp * ndofs)?; 1329 /// v.set_value(0.0); 1330 /// 1331 /// // Restrictions 1332 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1333 /// for i in 0..ne { 1334 /// indx[2 * i + 0] = i as i32; 1335 /// indx[2 * i + 1] = (i + 1) as i32; 1336 /// } 1337 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1338 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1339 /// for i in 0..ne { 1340 /// indu[p * i + 0] = (2 * i) as i32; 1341 /// indu[p * i + 1] = (2 * i + 1) as i32; 1342 /// indu[p * i + 2] = (2 * i + 2) as i32; 1343 /// } 1344 /// let ru = ceed.elem_restriction(ne, p, ncomp, ndofs, ncomp * ndofs, MemType::Host, &indu)?; 1345 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1346 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1347 /// 1348 /// // Bases 1349 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1350 /// let bu = ceed.basis_tensor_H1_Lagrange(1, ncomp, p, q, QuadMode::Gauss)?; 1351 /// 1352 /// // Build quadrature data 1353 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1354 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1355 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1356 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1357 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1358 /// .apply(&x, &mut qdata)?; 1359 /// 1360 /// // Mass operator 1361 /// let mut mass_2_comp = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1362 /// // Number of quadrature points 1363 /// let q = qdata.len(); 1364 /// 1365 /// // Iterate over quadrature points 1366 /// for i in 0..q { 1367 /// v[i + 0 * q] = u[i + 1 * q] * qdata[i]; 1368 /// v[i + 1 * q] = u[i + 0 * q] * qdata[i]; 1369 /// } 1370 /// 1371 /// // Return clean error code 1372 /// 0 1373 /// }; 1374 /// 1375 /// let qf_mass = ceed 1376 /// .q_function_interior(1, Box::new(mass_2_comp))? 1377 /// .input("u", 2, EvalMode::Interp)? 1378 /// .input("qdata", 1, EvalMode::None)? 1379 /// .output("v", 2, EvalMode::Interp)?; 1380 /// 1381 /// let op_mass = ceed 1382 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1383 /// .field("u", &ru, &bu, VectorOpt::Active)? 1384 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1385 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1386 /// 1387 /// // Diagonal 1388 /// let mut diag = ceed.vector(ncomp * ncomp * ndofs)?; 1389 /// diag.set_value(0.0); 1390 /// op_mass.linear_assemble_point_block_diagonal(&mut diag)?; 1391 /// 1392 /// // Manual diagonal computation 1393 /// let mut true_diag = ceed.vector(ncomp * ncomp * ndofs)?; 1394 /// true_diag.set_value(0.0)?; 1395 /// for i in 0..ndofs { 1396 /// for j in 0..ncomp { 1397 /// u.set_value(0.0); 1398 /// { 1399 /// let mut u_array = u.view_mut()?; 1400 /// u_array[i + j * ndofs] = 1.; 1401 /// } 1402 /// 1403 /// op_mass.apply(&u, &mut v)?; 1404 /// 1405 /// { 1406 /// let v_array = v.view()?; 1407 /// let mut true_array = true_diag.view_mut()?; 1408 /// for k in 0..ncomp { 1409 /// true_array[i * ncomp * ncomp + k * ncomp + j] = v_array[i + k * ndofs]; 1410 /// } 1411 /// } 1412 /// } 1413 /// } 1414 /// 1415 /// // Check 1416 /// diag.view()? 1417 /// .iter() 1418 /// .zip(true_diag.view()?.iter()) 1419 /// .for_each(|(computed, actual)| { 1420 /// assert!( 1421 /// (*computed - *actual).abs() < 10.0 * libceed::EPSILON, 1422 /// "Diagonal entry incorrect" 1423 /// ); 1424 /// }); 1425 /// # Ok(()) 1426 /// # } 1427 /// ``` 1428 pub fn linear_assemble_point_block_diagonal( 1429 &self, 1430 assembled: &mut Vector, 1431 ) -> crate::Result<i32> { 1432 self.op_core.linear_assemble_point_block_diagonal(assembled) 1433 } 1434 1435 /// Assemble the point block diagonal of a square linear Operator 1436 /// 1437 /// This sums into a Vector with the point block diagonal of a linear 1438 /// Operator. 1439 /// 1440 /// Note: Currently only non-composite Operators with a single field and 1441 /// composite Operators with single field sub-operators are supported. 1442 /// 1443 /// * `op` - Operator to assemble QFunction 1444 /// * `assembled` - Vector to store assembled CeedOperator point block 1445 /// diagonal, provided in row-major form with an 1446 /// `ncomp * ncomp` block at each node. The dimensions of 1447 /// this vector are derived from the active vector for 1448 /// the CeedOperator. The array has shape 1449 /// `[nodes, component out, component in]`. 1450 /// 1451 /// ``` 1452 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, EvalMode, MemType, QFunctionInputs, QFunctionOpt, QFunctionOutputs, QuadMode, Scalar, VectorOpt}; 1453 /// # fn main() -> libceed::Result<()> { 1454 /// # let ceed = libceed::Ceed::default_init(); 1455 /// let ne = 4; 1456 /// let p = 3; 1457 /// let q = 4; 1458 /// let ncomp = 2; 1459 /// let ndofs = p * ne - ne + 1; 1460 /// 1461 /// // Vectors 1462 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 1463 /// let mut qdata = ceed.vector(ne * q)?; 1464 /// qdata.set_value(0.0); 1465 /// let mut u = ceed.vector(ncomp * ndofs)?; 1466 /// u.set_value(1.0); 1467 /// let mut v = ceed.vector(ncomp * ndofs)?; 1468 /// v.set_value(0.0); 1469 /// 1470 /// // Restrictions 1471 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1472 /// for i in 0..ne { 1473 /// indx[2 * i + 0] = i as i32; 1474 /// indx[2 * i + 1] = (i + 1) as i32; 1475 /// } 1476 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1477 /// let mut indu: Vec<i32> = vec![0; p * ne]; 1478 /// for i in 0..ne { 1479 /// indu[p * i + 0] = (2 * i) as i32; 1480 /// indu[p * i + 1] = (2 * i + 1) as i32; 1481 /// indu[p * i + 2] = (2 * i + 2) as i32; 1482 /// } 1483 /// let ru = ceed.elem_restriction(ne, p, ncomp, ndofs, ncomp * ndofs, MemType::Host, &indu)?; 1484 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1485 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1486 /// 1487 /// // Bases 1488 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1489 /// let bu = ceed.basis_tensor_H1_Lagrange(1, ncomp, p, q, QuadMode::Gauss)?; 1490 /// 1491 /// // Build quadrature data 1492 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1493 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1494 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1495 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1496 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1497 /// .apply(&x, &mut qdata)?; 1498 /// 1499 /// // Mass operator 1500 /// let mut mass_2_comp = |[u, qdata, ..]: QFunctionInputs, [v, ..]: QFunctionOutputs| { 1501 /// // Number of quadrature points 1502 /// let q = qdata.len(); 1503 /// 1504 /// // Iterate over quadrature points 1505 /// for i in 0..q { 1506 /// v[i + 0 * q] = u[i + 1 * q] * qdata[i]; 1507 /// v[i + 1 * q] = u[i + 0 * q] * qdata[i]; 1508 /// } 1509 /// 1510 /// // Return clean error code 1511 /// 0 1512 /// }; 1513 /// 1514 /// let qf_mass = ceed 1515 /// .q_function_interior(1, Box::new(mass_2_comp))? 1516 /// .input("u", 2, EvalMode::Interp)? 1517 /// .input("qdata", 1, EvalMode::None)? 1518 /// .output("v", 2, EvalMode::Interp)?; 1519 /// 1520 /// let op_mass = ceed 1521 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1522 /// .field("u", &ru, &bu, VectorOpt::Active)? 1523 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1524 /// .field("v", &ru, &bu, VectorOpt::Active)?; 1525 /// 1526 /// // Diagonal 1527 /// let mut diag = ceed.vector(ncomp * ncomp * ndofs)?; 1528 /// diag.set_value(1.0); 1529 /// op_mass.linear_assemble_add_point_block_diagonal(&mut diag)?; 1530 /// 1531 /// // Manual diagonal computation 1532 /// let mut true_diag = ceed.vector(ncomp * ncomp * ndofs)?; 1533 /// true_diag.set_value(0.0)?; 1534 /// for i in 0..ndofs { 1535 /// for j in 0..ncomp { 1536 /// u.set_value(0.0); 1537 /// { 1538 /// let mut u_array = u.view_mut()?; 1539 /// u_array[i + j * ndofs] = 1.; 1540 /// } 1541 /// 1542 /// op_mass.apply(&u, &mut v)?; 1543 /// 1544 /// { 1545 /// let v_array = v.view()?; 1546 /// let mut true_array = true_diag.view_mut()?; 1547 /// for k in 0..ncomp { 1548 /// true_array[i * ncomp * ncomp + k * ncomp + j] = v_array[i + k * ndofs]; 1549 /// } 1550 /// } 1551 /// } 1552 /// } 1553 /// 1554 /// // Check 1555 /// diag.view()? 1556 /// .iter() 1557 /// .zip(true_diag.view()?.iter()) 1558 /// .for_each(|(computed, actual)| { 1559 /// assert!( 1560 /// (*computed - 1.0 - *actual).abs() < 10.0 * libceed::EPSILON, 1561 /// "Diagonal entry incorrect" 1562 /// ); 1563 /// }); 1564 /// # Ok(()) 1565 /// # } 1566 /// ``` 1567 pub fn linear_assemble_add_point_block_diagonal( 1568 &self, 1569 assembled: &mut Vector, 1570 ) -> crate::Result<i32> { 1571 self.op_core 1572 .linear_assemble_add_point_block_diagonal(assembled) 1573 } 1574 1575 /// Create a multigrid coarse Operator and level transfer Operators for a 1576 /// given Operator, creating the prolongation basis from the fine and 1577 /// coarse grid interpolation 1578 /// 1579 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1580 /// * `rstr_coarse` - Coarse grid restriction 1581 /// * `basis_coarse` - Coarse grid active vector basis 1582 /// 1583 /// ``` 1584 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 1585 /// # fn main() -> libceed::Result<()> { 1586 /// # let ceed = libceed::Ceed::default_init(); 1587 /// let ne = 15; 1588 /// let p_coarse = 3; 1589 /// let p_fine = 5; 1590 /// let q = 6; 1591 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1592 /// let ndofs_fine = p_fine * ne - ne + 1; 1593 /// 1594 /// // Vectors 1595 /// let x_array = (0..ne + 1) 1596 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1597 /// .collect::<Vec<Scalar>>(); 1598 /// let x = ceed.vector_from_slice(&x_array)?; 1599 /// let mut qdata = ceed.vector(ne * q)?; 1600 /// qdata.set_value(0.0); 1601 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1602 /// u_coarse.set_value(1.0); 1603 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1604 /// u_fine.set_value(1.0); 1605 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1606 /// v_coarse.set_value(0.0); 1607 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1608 /// v_fine.set_value(0.0); 1609 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1610 /// multiplicity.set_value(1.0); 1611 /// 1612 /// // Restrictions 1613 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1614 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1615 /// 1616 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1617 /// for i in 0..ne { 1618 /// indx[2 * i + 0] = i as i32; 1619 /// indx[2 * i + 1] = (i + 1) as i32; 1620 /// } 1621 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1622 /// 1623 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1624 /// for i in 0..ne { 1625 /// for j in 0..p_coarse { 1626 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1627 /// } 1628 /// } 1629 /// let ru_coarse = ceed.elem_restriction( 1630 /// ne, 1631 /// p_coarse, 1632 /// 1, 1633 /// 1, 1634 /// ndofs_coarse, 1635 /// MemType::Host, 1636 /// &indu_coarse, 1637 /// )?; 1638 /// 1639 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1640 /// for i in 0..ne { 1641 /// for j in 0..p_fine { 1642 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1643 /// } 1644 /// } 1645 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1646 /// 1647 /// // Bases 1648 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1649 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1650 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 1651 /// 1652 /// // Build quadrature data 1653 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1654 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1655 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1656 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1657 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1658 /// .apply(&x, &mut qdata)?; 1659 /// 1660 /// // Mass operator 1661 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1662 /// let op_mass_fine = ceed 1663 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1664 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 1665 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1666 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 1667 /// 1668 /// // Multigrid setup 1669 /// let (op_mass_coarse, op_prolong, op_restrict) = 1670 /// op_mass_fine.create_multigrid_level(&multiplicity, &ru_coarse, &bu_coarse)?; 1671 /// 1672 /// // Coarse problem 1673 /// u_coarse.set_value(1.0); 1674 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 1675 /// 1676 /// // Check 1677 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1678 /// assert!( 1679 /// (sum - 2.0).abs() < 50.0 * libceed::EPSILON, 1680 /// "Incorrect interval length computed" 1681 /// ); 1682 /// 1683 /// // Prolong 1684 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 1685 /// 1686 /// // Fine problem 1687 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 1688 /// 1689 /// // Check 1690 /// let sum: Scalar = v_fine.view()?.iter().sum(); 1691 /// assert!( 1692 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1693 /// "Incorrect interval length computed" 1694 /// ); 1695 /// 1696 /// // Restrict 1697 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 1698 /// 1699 /// // Check 1700 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1701 /// assert!( 1702 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1703 /// "Incorrect interval length computed" 1704 /// ); 1705 /// # Ok(()) 1706 /// # } 1707 /// ``` 1708 pub fn create_multigrid_level<'b>( 1709 &self, 1710 p_mult_fine: &Vector, 1711 rstr_coarse: &ElemRestriction, 1712 basis_coarse: &Basis, 1713 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 1714 let mut ptr_coarse = std::ptr::null_mut(); 1715 let mut ptr_prolong = std::ptr::null_mut(); 1716 let mut ptr_restrict = std::ptr::null_mut(); 1717 self.op_core.check_error(unsafe { 1718 bind_ceed::CeedOperatorMultigridLevelCreate( 1719 self.op_core.ptr, 1720 p_mult_fine.ptr, 1721 rstr_coarse.ptr, 1722 basis_coarse.ptr, 1723 &mut ptr_coarse, 1724 &mut ptr_prolong, 1725 &mut ptr_restrict, 1726 ) 1727 })?; 1728 let op_coarse = unsafe { Operator::from_raw(ptr_coarse)? }; 1729 let op_prolong = unsafe { Operator::from_raw(ptr_prolong)? }; 1730 let op_restrict = unsafe { Operator::from_raw(ptr_restrict)? }; 1731 Ok((op_coarse, op_prolong, op_restrict)) 1732 } 1733 1734 /// Create a multigrid coarse Operator and level transfer Operators for a 1735 /// given Operator with a tensor basis for the active basis 1736 /// 1737 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1738 /// * `rstr_coarse` - Coarse grid restriction 1739 /// * `basis_coarse` - Coarse grid active vector basis 1740 /// * `interp_c_to_f` - Matrix for coarse to fine 1741 /// 1742 /// ``` 1743 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, EvalMode, MemType, QFunctionOpt, QuadMode, Scalar, TransposeMode, VectorOpt}; 1744 /// # fn main() -> libceed::Result<()> { 1745 /// # let ceed = libceed::Ceed::default_init(); 1746 /// let ne = 15; 1747 /// let p_coarse = 3; 1748 /// let p_fine = 5; 1749 /// let q = 6; 1750 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1751 /// let ndofs_fine = p_fine * ne - ne + 1; 1752 /// 1753 /// // Vectors 1754 /// let x_array = (0..ne + 1) 1755 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1756 /// .collect::<Vec<Scalar>>(); 1757 /// let x = ceed.vector_from_slice(&x_array)?; 1758 /// let mut qdata = ceed.vector(ne * q)?; 1759 /// qdata.set_value(0.0); 1760 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1761 /// u_coarse.set_value(1.0); 1762 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1763 /// u_fine.set_value(1.0); 1764 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1765 /// v_coarse.set_value(0.0); 1766 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1767 /// v_fine.set_value(0.0); 1768 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1769 /// multiplicity.set_value(1.0); 1770 /// 1771 /// // Restrictions 1772 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1773 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1774 /// 1775 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1776 /// for i in 0..ne { 1777 /// indx[2 * i + 0] = i as i32; 1778 /// indx[2 * i + 1] = (i + 1) as i32; 1779 /// } 1780 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1781 /// 1782 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1783 /// for i in 0..ne { 1784 /// for j in 0..p_coarse { 1785 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1786 /// } 1787 /// } 1788 /// let ru_coarse = ceed.elem_restriction( 1789 /// ne, 1790 /// p_coarse, 1791 /// 1, 1792 /// 1, 1793 /// ndofs_coarse, 1794 /// MemType::Host, 1795 /// &indu_coarse, 1796 /// )?; 1797 /// 1798 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1799 /// for i in 0..ne { 1800 /// for j in 0..p_fine { 1801 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1802 /// } 1803 /// } 1804 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1805 /// 1806 /// // Bases 1807 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1808 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1809 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 1810 /// 1811 /// // Build quadrature data 1812 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 1813 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 1814 /// .field("dx", &rx, &bx, VectorOpt::Active)? 1815 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 1816 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 1817 /// .apply(&x, &mut qdata)?; 1818 /// 1819 /// // Mass operator 1820 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 1821 /// let op_mass_fine = ceed 1822 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 1823 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 1824 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 1825 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 1826 /// 1827 /// // Multigrid setup 1828 /// let mut interp_c_to_f: Vec<Scalar> = vec![0.; p_coarse * p_fine]; 1829 /// { 1830 /// let mut coarse = ceed.vector(p_coarse)?; 1831 /// let mut fine = ceed.vector(p_fine)?; 1832 /// let basis_c_to_f = 1833 /// ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, p_fine, QuadMode::GaussLobatto)?; 1834 /// for i in 0..p_coarse { 1835 /// coarse.set_value(0.0); 1836 /// { 1837 /// let mut array = coarse.view_mut()?; 1838 /// array[i] = 1.; 1839 /// } 1840 /// basis_c_to_f.apply( 1841 /// 1, 1842 /// TransposeMode::NoTranspose, 1843 /// EvalMode::Interp, 1844 /// &coarse, 1845 /// &mut fine, 1846 /// )?; 1847 /// let array = fine.view()?; 1848 /// for j in 0..p_fine { 1849 /// interp_c_to_f[j * p_coarse + i] = array[j]; 1850 /// } 1851 /// } 1852 /// } 1853 /// let (op_mass_coarse, op_prolong, op_restrict) = op_mass_fine.create_multigrid_level_tensor_H1( 1854 /// &multiplicity, 1855 /// &ru_coarse, 1856 /// &bu_coarse, 1857 /// &interp_c_to_f, 1858 /// )?; 1859 /// 1860 /// // Coarse problem 1861 /// u_coarse.set_value(1.0); 1862 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 1863 /// 1864 /// // Check 1865 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1866 /// assert!( 1867 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 1868 /// "Incorrect interval length computed" 1869 /// ); 1870 /// 1871 /// // Prolong 1872 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 1873 /// 1874 /// // Fine problem 1875 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 1876 /// 1877 /// // Check 1878 /// let sum: Scalar = v_fine.view()?.iter().sum(); 1879 /// assert!( 1880 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1881 /// "Incorrect interval length computed" 1882 /// ); 1883 /// 1884 /// // Restrict 1885 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 1886 /// 1887 /// // Check 1888 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 1889 /// assert!( 1890 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 1891 /// "Incorrect interval length computed" 1892 /// ); 1893 /// # Ok(()) 1894 /// # } 1895 /// ``` 1896 pub fn create_multigrid_level_tensor_H1<'b>( 1897 &self, 1898 p_mult_fine: &Vector, 1899 rstr_coarse: &ElemRestriction, 1900 basis_coarse: &Basis, 1901 interpCtoF: &[crate::Scalar], 1902 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 1903 let mut ptr_coarse = std::ptr::null_mut(); 1904 let mut ptr_prolong = std::ptr::null_mut(); 1905 let mut ptr_restrict = std::ptr::null_mut(); 1906 self.op_core.check_error(unsafe { 1907 bind_ceed::CeedOperatorMultigridLevelCreateTensorH1( 1908 self.op_core.ptr, 1909 p_mult_fine.ptr, 1910 rstr_coarse.ptr, 1911 basis_coarse.ptr, 1912 interpCtoF.as_ptr(), 1913 &mut ptr_coarse, 1914 &mut ptr_prolong, 1915 &mut ptr_restrict, 1916 ) 1917 })?; 1918 let op_coarse = unsafe { Operator::from_raw(ptr_coarse)? }; 1919 let op_prolong = unsafe { Operator::from_raw(ptr_prolong)? }; 1920 let op_restrict = unsafe { Operator::from_raw(ptr_restrict)? }; 1921 Ok((op_coarse, op_prolong, op_restrict)) 1922 } 1923 1924 /// Create a multigrid coarse Operator and level transfer Operators for a 1925 /// given Operator with a non-tensor basis for the active basis 1926 /// 1927 /// * `p_mult_fine` - Lvector multiplicity in parallel gather/scatter 1928 /// * `rstr_coarse` - Coarse grid restriction 1929 /// * `basis_coarse` - Coarse grid active vector basis 1930 /// * `interp_c_to_f` - Matrix for coarse to fine 1931 /// 1932 /// ``` 1933 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, EvalMode, MemType, QFunctionOpt, QuadMode, Scalar, TransposeMode, VectorOpt}; 1934 /// # fn main() -> libceed::Result<()> { 1935 /// # let ceed = libceed::Ceed::default_init(); 1936 /// let ne = 15; 1937 /// let p_coarse = 3; 1938 /// let p_fine = 5; 1939 /// let q = 6; 1940 /// let ndofs_coarse = p_coarse * ne - ne + 1; 1941 /// let ndofs_fine = p_fine * ne - ne + 1; 1942 /// 1943 /// // Vectors 1944 /// let x_array = (0..ne + 1) 1945 /// .map(|i| 2.0 * i as Scalar / ne as Scalar - 1.0) 1946 /// .collect::<Vec<Scalar>>(); 1947 /// let x = ceed.vector_from_slice(&x_array)?; 1948 /// let mut qdata = ceed.vector(ne * q)?; 1949 /// qdata.set_value(0.0); 1950 /// let mut u_coarse = ceed.vector(ndofs_coarse)?; 1951 /// u_coarse.set_value(1.0); 1952 /// let mut u_fine = ceed.vector(ndofs_fine)?; 1953 /// u_fine.set_value(1.0); 1954 /// let mut v_coarse = ceed.vector(ndofs_coarse)?; 1955 /// v_coarse.set_value(0.0); 1956 /// let mut v_fine = ceed.vector(ndofs_fine)?; 1957 /// v_fine.set_value(0.0); 1958 /// let mut multiplicity = ceed.vector(ndofs_fine)?; 1959 /// multiplicity.set_value(1.0); 1960 /// 1961 /// // Restrictions 1962 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 1963 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 1964 /// 1965 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 1966 /// for i in 0..ne { 1967 /// indx[2 * i + 0] = i as i32; 1968 /// indx[2 * i + 1] = (i + 1) as i32; 1969 /// } 1970 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 1971 /// 1972 /// let mut indu_coarse: Vec<i32> = vec![0; p_coarse * ne]; 1973 /// for i in 0..ne { 1974 /// for j in 0..p_coarse { 1975 /// indu_coarse[p_coarse * i + j] = (i + j) as i32; 1976 /// } 1977 /// } 1978 /// let ru_coarse = ceed.elem_restriction( 1979 /// ne, 1980 /// p_coarse, 1981 /// 1, 1982 /// 1, 1983 /// ndofs_coarse, 1984 /// MemType::Host, 1985 /// &indu_coarse, 1986 /// )?; 1987 /// 1988 /// let mut indu_fine: Vec<i32> = vec![0; p_fine * ne]; 1989 /// for i in 0..ne { 1990 /// for j in 0..p_fine { 1991 /// indu_fine[p_fine * i + j] = (i + j) as i32; 1992 /// } 1993 /// } 1994 /// let ru_fine = ceed.elem_restriction(ne, p_fine, 1, 1, ndofs_fine, MemType::Host, &indu_fine)?; 1995 /// 1996 /// // Bases 1997 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 1998 /// let bu_coarse = ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, q, QuadMode::Gauss)?; 1999 /// let bu_fine = ceed.basis_tensor_H1_Lagrange(1, 1, p_fine, q, QuadMode::Gauss)?; 2000 /// 2001 /// // Build quadrature data 2002 /// let qf_build = ceed.q_function_interior_by_name("Mass1DBuild")?; 2003 /// ceed.operator(&qf_build, QFunctionOpt::None, QFunctionOpt::None)? 2004 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2005 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2006 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2007 /// .apply(&x, &mut qdata)?; 2008 /// 2009 /// // Mass operator 2010 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2011 /// let op_mass_fine = ceed 2012 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2013 /// .field("u", &ru_fine, &bu_fine, VectorOpt::Active)? 2014 /// .field("qdata", &rq, BasisOpt::None, &qdata)? 2015 /// .field("v", &ru_fine, &bu_fine, VectorOpt::Active)?; 2016 /// 2017 /// // Multigrid setup 2018 /// let mut interp_c_to_f: Vec<Scalar> = vec![0.; p_coarse * p_fine]; 2019 /// { 2020 /// let mut coarse = ceed.vector(p_coarse)?; 2021 /// let mut fine = ceed.vector(p_fine)?; 2022 /// let basis_c_to_f = 2023 /// ceed.basis_tensor_H1_Lagrange(1, 1, p_coarse, p_fine, QuadMode::GaussLobatto)?; 2024 /// for i in 0..p_coarse { 2025 /// coarse.set_value(0.0); 2026 /// { 2027 /// let mut array = coarse.view_mut()?; 2028 /// array[i] = 1.; 2029 /// } 2030 /// basis_c_to_f.apply( 2031 /// 1, 2032 /// TransposeMode::NoTranspose, 2033 /// EvalMode::Interp, 2034 /// &coarse, 2035 /// &mut fine, 2036 /// )?; 2037 /// let array = fine.view()?; 2038 /// for j in 0..p_fine { 2039 /// interp_c_to_f[j * p_coarse + i] = array[j]; 2040 /// } 2041 /// } 2042 /// } 2043 /// let (op_mass_coarse, op_prolong, op_restrict) = op_mass_fine.create_multigrid_level_H1( 2044 /// &multiplicity, 2045 /// &ru_coarse, 2046 /// &bu_coarse, 2047 /// &interp_c_to_f, 2048 /// )?; 2049 /// 2050 /// // Coarse problem 2051 /// u_coarse.set_value(1.0); 2052 /// op_mass_coarse.apply(&u_coarse, &mut v_coarse)?; 2053 /// 2054 /// // Check 2055 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 2056 /// assert!( 2057 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 2058 /// "Incorrect interval length computed" 2059 /// ); 2060 /// 2061 /// // Prolong 2062 /// op_prolong.apply(&u_coarse, &mut u_fine)?; 2063 /// 2064 /// // Fine problem 2065 /// op_mass_fine.apply(&u_fine, &mut v_fine)?; 2066 /// 2067 /// // Check 2068 /// let sum: Scalar = v_fine.view()?.iter().sum(); 2069 /// assert!( 2070 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 2071 /// "Incorrect interval length computed" 2072 /// ); 2073 /// 2074 /// // Restrict 2075 /// op_restrict.apply(&v_fine, &mut v_coarse)?; 2076 /// 2077 /// // Check 2078 /// let sum: Scalar = v_coarse.view()?.iter().sum(); 2079 /// assert!( 2080 /// (sum - 2.0).abs() < 200.0 * libceed::EPSILON, 2081 /// "Incorrect interval length computed" 2082 /// ); 2083 /// # Ok(()) 2084 /// # } 2085 /// ``` 2086 pub fn create_multigrid_level_H1<'b>( 2087 &self, 2088 p_mult_fine: &Vector, 2089 rstr_coarse: &ElemRestriction, 2090 basis_coarse: &Basis, 2091 interpCtoF: &[crate::Scalar], 2092 ) -> crate::Result<(Operator<'b>, Operator<'b>, Operator<'b>)> { 2093 let mut ptr_coarse = std::ptr::null_mut(); 2094 let mut ptr_prolong = std::ptr::null_mut(); 2095 let mut ptr_restrict = std::ptr::null_mut(); 2096 self.op_core.check_error(unsafe { 2097 bind_ceed::CeedOperatorMultigridLevelCreateH1( 2098 self.op_core.ptr, 2099 p_mult_fine.ptr, 2100 rstr_coarse.ptr, 2101 basis_coarse.ptr, 2102 interpCtoF.as_ptr(), 2103 &mut ptr_coarse, 2104 &mut ptr_prolong, 2105 &mut ptr_restrict, 2106 ) 2107 })?; 2108 let op_coarse = unsafe { Operator::from_raw(ptr_coarse)? }; 2109 let op_prolong = unsafe { Operator::from_raw(ptr_prolong)? }; 2110 let op_restrict = unsafe { Operator::from_raw(ptr_restrict)? }; 2111 Ok((op_coarse, op_prolong, op_restrict)) 2112 } 2113 } 2114 2115 // ----------------------------------------------------------------------------- 2116 // Composite Operator 2117 // ----------------------------------------------------------------------------- 2118 impl<'a> CompositeOperator<'a> { 2119 // Constructor 2120 pub fn create(ceed: &crate::Ceed) -> crate::Result<Self> { 2121 let mut ptr = std::ptr::null_mut(); 2122 ceed.check_error(unsafe { bind_ceed::CeedCompositeOperatorCreate(ceed.ptr, &mut ptr) })?; 2123 Ok(Self { 2124 op_core: OperatorCore { 2125 ptr, 2126 _lifeline: PhantomData, 2127 }, 2128 }) 2129 } 2130 2131 /// Set name for CompositeOperator printing 2132 /// 2133 /// * 'name' - Name to set 2134 /// 2135 /// ``` 2136 /// # use libceed::{prelude::*, BasisOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 2137 /// # fn main() -> libceed::Result<()> { 2138 /// # let ceed = libceed::Ceed::default_init(); 2139 /// 2140 /// // Sub operator field arguments 2141 /// let ne = 3; 2142 /// let q = 4_usize; 2143 /// let mut ind: Vec<i32> = vec![0; 2 * ne]; 2144 /// for i in 0..ne { 2145 /// ind[2 * i + 0] = i as i32; 2146 /// ind[2 * i + 1] = (i + 1) as i32; 2147 /// } 2148 /// let r = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &ind)?; 2149 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2150 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2151 /// 2152 /// let b = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2153 /// 2154 /// let qdata_mass = ceed.vector(q * ne)?; 2155 /// let qdata_diff = ceed.vector(q * ne)?; 2156 /// 2157 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2158 /// let op_mass = ceed 2159 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2160 /// .name("Mass term")? 2161 /// .field("u", &r, &b, VectorOpt::Active)? 2162 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2163 /// .field("v", &r, &b, VectorOpt::Active)?; 2164 /// 2165 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2166 /// let op_diff = ceed 2167 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2168 /// .name("Poisson term")? 2169 /// .field("du", &r, &b, VectorOpt::Active)? 2170 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2171 /// .field("dv", &r, &b, VectorOpt::Active)?; 2172 /// 2173 /// let op = ceed 2174 /// .composite_operator()? 2175 /// .name("Screened Poisson")? 2176 /// .sub_operator(&op_mass)? 2177 /// .sub_operator(&op_diff)?; 2178 /// # Ok(()) 2179 /// # } 2180 /// ``` 2181 #[allow(unused_mut)] 2182 pub fn name(mut self, name: &str) -> crate::Result<Self> { 2183 self.op_core.name(name)?; 2184 Ok(self) 2185 } 2186 2187 /// Apply Operator to a vector 2188 /// 2189 /// * `input` - Inpuht Vector 2190 /// * `output` - Output Vector 2191 /// 2192 /// ``` 2193 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 2194 /// # fn main() -> libceed::Result<()> { 2195 /// # let ceed = libceed::Ceed::default_init(); 2196 /// let ne = 4; 2197 /// let p = 3; 2198 /// let q = 4; 2199 /// let ndofs = p * ne - ne + 1; 2200 /// 2201 /// // Vectors 2202 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 2203 /// let mut qdata_mass = ceed.vector(ne * q)?; 2204 /// qdata_mass.set_value(0.0); 2205 /// let mut qdata_diff = ceed.vector(ne * q)?; 2206 /// qdata_diff.set_value(0.0); 2207 /// let mut u = ceed.vector(ndofs)?; 2208 /// u.set_value(1.0); 2209 /// let mut v = ceed.vector(ndofs)?; 2210 /// v.set_value(0.0); 2211 /// 2212 /// // Restrictions 2213 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2214 /// for i in 0..ne { 2215 /// indx[2 * i + 0] = i as i32; 2216 /// indx[2 * i + 1] = (i + 1) as i32; 2217 /// } 2218 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2219 /// let mut indu: Vec<i32> = vec![0; p * ne]; 2220 /// for i in 0..ne { 2221 /// indu[p * i + 0] = i as i32; 2222 /// indu[p * i + 1] = (i + 1) as i32; 2223 /// indu[p * i + 2] = (i + 2) as i32; 2224 /// } 2225 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 2226 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2227 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2228 /// 2229 /// // Bases 2230 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2231 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 2232 /// 2233 /// // Build quadrature data 2234 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2235 /// ceed.operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2236 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2237 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2238 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2239 /// .apply(&x, &mut qdata_mass)?; 2240 /// 2241 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2242 /// ceed.operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2243 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2244 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2245 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2246 /// .apply(&x, &mut qdata_diff)?; 2247 /// 2248 /// // Application operator 2249 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2250 /// let op_mass = ceed 2251 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2252 /// .field("u", &ru, &bu, VectorOpt::Active)? 2253 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2254 /// .field("v", &ru, &bu, VectorOpt::Active)?; 2255 /// 2256 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2257 /// let op_diff = ceed 2258 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2259 /// .field("du", &ru, &bu, VectorOpt::Active)? 2260 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2261 /// .field("dv", &ru, &bu, VectorOpt::Active)?; 2262 /// 2263 /// let op_composite = ceed 2264 /// .composite_operator()? 2265 /// .sub_operator(&op_mass)? 2266 /// .sub_operator(&op_diff)?; 2267 /// 2268 /// v.set_value(0.0); 2269 /// op_composite.apply(&u, &mut v)?; 2270 /// 2271 /// // Check 2272 /// let sum: Scalar = v.view()?.iter().sum(); 2273 /// assert!( 2274 /// (sum - 2.0).abs() < 10.0 * libceed::EPSILON, 2275 /// "Incorrect interval length computed" 2276 /// ); 2277 /// # Ok(()) 2278 /// # } 2279 /// ``` 2280 pub fn apply(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 2281 self.op_core.apply(input, output) 2282 } 2283 2284 /// Apply Operator to a vector and add result to output vector 2285 /// 2286 /// * `input` - Input Vector 2287 /// * `output` - Output Vector 2288 /// 2289 /// ``` 2290 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, Scalar, VectorOpt}; 2291 /// # fn main() -> libceed::Result<()> { 2292 /// # let ceed = libceed::Ceed::default_init(); 2293 /// let ne = 4; 2294 /// let p = 3; 2295 /// let q = 4; 2296 /// let ndofs = p * ne - ne + 1; 2297 /// 2298 /// // Vectors 2299 /// let x = ceed.vector_from_slice(&[-1., -0.5, 0.0, 0.5, 1.0])?; 2300 /// let mut qdata_mass = ceed.vector(ne * q)?; 2301 /// qdata_mass.set_value(0.0); 2302 /// let mut qdata_diff = ceed.vector(ne * q)?; 2303 /// qdata_diff.set_value(0.0); 2304 /// let mut u = ceed.vector(ndofs)?; 2305 /// u.set_value(1.0); 2306 /// let mut v = ceed.vector(ndofs)?; 2307 /// v.set_value(0.0); 2308 /// 2309 /// // Restrictions 2310 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2311 /// for i in 0..ne { 2312 /// indx[2 * i + 0] = i as i32; 2313 /// indx[2 * i + 1] = (i + 1) as i32; 2314 /// } 2315 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2316 /// let mut indu: Vec<i32> = vec![0; p * ne]; 2317 /// for i in 0..ne { 2318 /// indu[p * i + 0] = i as i32; 2319 /// indu[p * i + 1] = (i + 1) as i32; 2320 /// indu[p * i + 2] = (i + 2) as i32; 2321 /// } 2322 /// let ru = ceed.elem_restriction(ne, 3, 1, 1, ndofs, MemType::Host, &indu)?; 2323 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2324 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2325 /// 2326 /// // Bases 2327 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2328 /// let bu = ceed.basis_tensor_H1_Lagrange(1, 1, p, q, QuadMode::Gauss)?; 2329 /// 2330 /// // Build quadrature data 2331 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2332 /// ceed.operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2333 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2334 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2335 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2336 /// .apply(&x, &mut qdata_mass)?; 2337 /// 2338 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2339 /// ceed.operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2340 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2341 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2342 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)? 2343 /// .apply(&x, &mut qdata_diff)?; 2344 /// 2345 /// // Application operator 2346 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2347 /// let op_mass = ceed 2348 /// .operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)? 2349 /// .field("u", &ru, &bu, VectorOpt::Active)? 2350 /// .field("qdata", &rq, BasisOpt::None, &qdata_mass)? 2351 /// .field("v", &ru, &bu, VectorOpt::Active)?; 2352 /// 2353 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2354 /// let op_diff = ceed 2355 /// .operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)? 2356 /// .field("du", &ru, &bu, VectorOpt::Active)? 2357 /// .field("qdata", &rq, BasisOpt::None, &qdata_diff)? 2358 /// .field("dv", &ru, &bu, VectorOpt::Active)?; 2359 /// 2360 /// let op_composite = ceed 2361 /// .composite_operator()? 2362 /// .sub_operator(&op_mass)? 2363 /// .sub_operator(&op_diff)?; 2364 /// 2365 /// v.set_value(1.0); 2366 /// op_composite.apply_add(&u, &mut v)?; 2367 /// 2368 /// // Check 2369 /// let sum: Scalar = v.view()?.iter().sum(); 2370 /// assert!( 2371 /// (sum - (2.0 + ndofs as Scalar)).abs() < 10.0 * libceed::EPSILON, 2372 /// "Incorrect interval length computed" 2373 /// ); 2374 /// # Ok(()) 2375 /// # } 2376 /// ``` 2377 pub fn apply_add(&self, input: &Vector, output: &mut Vector) -> crate::Result<i32> { 2378 self.op_core.apply_add(input, output) 2379 } 2380 2381 /// Add a sub-Operator to a Composite Operator 2382 /// 2383 /// * `subop` - Sub-Operator 2384 /// 2385 /// ``` 2386 /// # use libceed::{prelude::*, QFunctionOpt}; 2387 /// # fn main() -> libceed::Result<()> { 2388 /// # let ceed = libceed::Ceed::default_init(); 2389 /// let mut op = ceed.composite_operator()?; 2390 /// 2391 /// let qf_mass = ceed.q_function_interior_by_name("MassApply")?; 2392 /// let op_mass = ceed.operator(&qf_mass, QFunctionOpt::None, QFunctionOpt::None)?; 2393 /// op = op.sub_operator(&op_mass)?; 2394 /// 2395 /// let qf_diff = ceed.q_function_interior_by_name("Poisson1DApply")?; 2396 /// let op_diff = ceed.operator(&qf_diff, QFunctionOpt::None, QFunctionOpt::None)?; 2397 /// op = op.sub_operator(&op_diff)?; 2398 /// # Ok(()) 2399 /// # } 2400 /// ``` 2401 #[allow(unused_mut)] 2402 pub fn sub_operator(mut self, subop: &Operator) -> crate::Result<Self> { 2403 self.op_core.check_error(unsafe { 2404 bind_ceed::CeedCompositeOperatorAddSub(self.op_core.ptr, subop.op_core.ptr) 2405 })?; 2406 Ok(self) 2407 } 2408 2409 /// Check if CompositeOperator is setup correctly 2410 /// 2411 /// ``` 2412 /// # use libceed::{prelude::*, BasisOpt, ElemRestrictionOpt, MemType, QFunctionOpt, QuadMode, VectorOpt}; 2413 /// # fn main() -> libceed::Result<()> { 2414 /// # let ceed = libceed::Ceed::default_init(); 2415 /// let ne = 4; 2416 /// let p = 3; 2417 /// let q = 4; 2418 /// let ndofs = p * ne - ne + 1; 2419 /// 2420 /// // Restrictions 2421 /// let mut indx: Vec<i32> = vec![0; 2 * ne]; 2422 /// for i in 0..ne { 2423 /// indx[2 * i + 0] = i as i32; 2424 /// indx[2 * i + 1] = (i + 1) as i32; 2425 /// } 2426 /// let rx = ceed.elem_restriction(ne, 2, 1, 1, ne + 1, MemType::Host, &indx)?; 2427 /// let strides: [i32; 3] = [1, q as i32, q as i32]; 2428 /// let rq = ceed.strided_elem_restriction(ne, q, 1, q * ne, strides)?; 2429 /// 2430 /// // Bases 2431 /// let bx = ceed.basis_tensor_H1_Lagrange(1, 1, 2, q, QuadMode::Gauss)?; 2432 /// 2433 /// // Build quadrature data 2434 /// let qf_build_mass = ceed.q_function_interior_by_name("Mass1DBuild")?; 2435 /// let op_build_mass = ceed 2436 /// .operator(&qf_build_mass, QFunctionOpt::None, QFunctionOpt::None)? 2437 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2438 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2439 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 2440 /// 2441 /// let qf_build_diff = ceed.q_function_interior_by_name("Poisson1DBuild")?; 2442 /// let op_build_diff = ceed 2443 /// .operator(&qf_build_diff, QFunctionOpt::None, QFunctionOpt::None)? 2444 /// .field("dx", &rx, &bx, VectorOpt::Active)? 2445 /// .field("weights", ElemRestrictionOpt::None, &bx, VectorOpt::None)? 2446 /// .field("qdata", &rq, BasisOpt::None, VectorOpt::Active)?; 2447 /// 2448 /// let op_build = ceed 2449 /// .composite_operator()? 2450 /// .sub_operator(&op_build_mass)? 2451 /// .sub_operator(&op_build_diff)?; 2452 /// 2453 /// // Check 2454 /// assert!(op_build.check().is_ok(), "Operator setup incorrectly"); 2455 /// # Ok(()) 2456 /// # } 2457 /// ``` 2458 pub fn check(self) -> crate::Result<Self> { 2459 self.op_core.check()?; 2460 Ok(self) 2461 } 2462 2463 /// Assemble the diagonal of a square linear Operator 2464 /// 2465 /// This overwrites a Vector with the diagonal of a linear Operator. 2466 /// 2467 /// Note: Currently only non-composite Operators with a single field and 2468 /// composite Operators with single field sub-operators are supported. 2469 /// 2470 /// * `op` - Operator to assemble QFunction 2471 /// * `assembled` - Vector to store assembled Operator diagonal 2472 pub fn linear_assemble_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 2473 self.op_core.linear_assemble_add_diagonal(assembled) 2474 } 2475 2476 /// Assemble the point block diagonal of a square linear Operator 2477 /// 2478 /// This overwrites a Vector with the point block diagonal of a linear 2479 /// Operator. 2480 /// 2481 /// Note: Currently only non-composite Operators with a single field and 2482 /// composite Operators with single field sub-operators are supported. 2483 /// 2484 /// * `op` - Operator to assemble QFunction 2485 /// * `assembled` - Vector to store assembled Operator diagonal 2486 pub fn linear_assemble_add_diagonal(&self, assembled: &mut Vector) -> crate::Result<i32> { 2487 self.op_core.linear_assemble_add_diagonal(assembled) 2488 } 2489 2490 /// Assemble the diagonal of a square linear Operator 2491 /// 2492 /// This overwrites a Vector with the diagonal of a linear Operator. 2493 /// 2494 /// Note: Currently only non-composite Operators with a single field and 2495 /// composite Operators with single field sub-operators are supported. 2496 /// 2497 /// * `op` - Operator to assemble QFunction 2498 /// * `assembled` - Vector to store assembled CeedOperator point block 2499 /// diagonal, provided in row-major form with an 2500 /// `ncomp * ncomp` block at each node. The dimensions of 2501 /// this vector are derived from the active vector for 2502 /// the CeedOperator. The array has shape 2503 /// `[nodes, component out, component in]`. 2504 pub fn linear_assemble_point_block_diagonal( 2505 &self, 2506 assembled: &mut Vector, 2507 ) -> crate::Result<i32> { 2508 self.op_core.linear_assemble_add_diagonal(assembled) 2509 } 2510 2511 /// Assemble the diagonal of a square linear Operator 2512 /// 2513 /// This sums into a Vector with the diagonal of a linear Operator. 2514 /// 2515 /// Note: Currently only non-composite Operators with a single field and 2516 /// composite Operators with single field sub-operators are supported. 2517 /// 2518 /// * `op` - Operator to assemble QFunction 2519 /// * `assembled` - Vector to store assembled CeedOperator point block 2520 /// diagonal, provided in row-major form with an 2521 /// `ncomp * ncomp` block at each node. The dimensions of 2522 /// this vector are derived from the active vector for 2523 /// the CeedOperator. The array has shape 2524 /// `[nodes, component out, component in]`. 2525 pub fn linear_assemble_add_point_block_diagonal( 2526 &self, 2527 assembled: &mut Vector, 2528 ) -> crate::Result<i32> { 2529 self.op_core.linear_assemble_add_diagonal(assembled) 2530 } 2531 } 2532 2533 // ----------------------------------------------------------------------------- 2534