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