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