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