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