1 // Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 use libceed::{Ceed, ElemRestriction, Vector}; 9 use std::convert::TryInto; 10 11 // ---------------------------------------------------------------------------- 12 // Determine problem size in each dimension from size and dimenison 13 // ---------------------------------------------------------------------------- 14 pub fn cartesian_mesh_size(dim: usize, solution_degree: usize, problem_size: i64) -> [usize; 3] { 15 // Use the approximate formula: 16 // prob_size ~ num_elem * degree^dim 17 let mut num_elem = problem_size / solution_degree.pow(dim as u32) as i64; 18 let mut s = 0; // find s: num_elem / 2 < 2^s <= num_elem 19 while num_elem > 1 { 20 num_elem /= 2; 21 s += 1; 22 } 23 24 // Size per dimension 25 let mut r = s % dim; 26 let xyz: [usize; 3] = (0..3) 27 .map(|_| -> usize { 28 let mut sd = s / dim; 29 if r > 0 { 30 sd += 1; 31 r -= 1; 32 } 33 1 << sd 34 }) 35 .collect::<Vec<usize>>() 36 .try_into() 37 .unwrap(); 38 xyz 39 } 40 41 // ---------------------------------------------------------------------------- 42 // Build element restriction objects for the mesh 43 // ---------------------------------------------------------------------------- 44 pub fn build_cartesian_restriction( 45 ceed: &Ceed, 46 dim: usize, 47 num_xyz: [usize; 3], 48 degree: usize, 49 num_comp: usize, 50 num_qpts: usize, 51 ) -> libceed::Result<(ElemRestriction, ElemRestriction)> { 52 let p = degree + 1; 53 let num_nodes = p.pow(dim as u32); // number of nodes per element 54 let elem_qpts = num_qpts.pow(dim as u32); // number of quadrature pts per element 55 56 // Problem dimensions 57 let mut num_d = [0; 3]; 58 let mut num_elem = 1; 59 let mut scalar_size = 1; 60 for d in 0..dim { 61 num_elem *= num_xyz[d]; 62 num_d[d] = num_xyz[d] * (p - 1) + 1; 63 scalar_size *= num_d[d]; 64 } 65 66 // elem: 0 1 n-1 67 // |---*-...-*---|---*-...-*---|- ... -|--...--| 68 // nodes: 0 1 p-1 p p+1 2*p n*p 69 let mut elem_nodes = vec![0; num_elem * num_nodes]; 70 for e in 0..num_elem { 71 let mut e_xyz = [1; 3]; 72 let mut re = e; 73 for d in 0..dim { 74 e_xyz[d] = re % num_xyz[d]; 75 re /= num_xyz[d]; 76 } 77 let loc_offset = e * num_nodes; 78 for loc_nodes in 0..num_nodes { 79 let mut global_nodes = 0; 80 let mut global_nodes_stride = 1; 81 let mut r_nodes = loc_nodes; 82 for d in 0..dim { 83 global_nodes += (e_xyz[d] * (p - 1) + r_nodes % p) * global_nodes_stride; 84 global_nodes_stride *= num_d[d]; 85 r_nodes /= p; 86 } 87 elem_nodes[loc_offset + loc_nodes] = global_nodes as i32; 88 } 89 } 90 91 // Mesh/solution data restriction 92 let rstr = ceed.elem_restriction( 93 num_elem, 94 num_nodes, 95 num_comp, 96 scalar_size, 97 num_comp * scalar_size, 98 libceed::MemType::Host, 99 &elem_nodes, 100 )?; 101 102 // Quadrature data restriction 103 let rstr_qdata = ceed.strided_elem_restriction( 104 num_elem, 105 elem_qpts, 106 num_comp, 107 num_comp * elem_qpts * num_elem, 108 libceed::CEED_STRIDES_BACKEND, 109 )?; 110 Ok((rstr, rstr_qdata)) 111 } 112 113 // ---------------------------------------------------------------------------- 114 // Set mesh coordinates 115 // ---------------------------------------------------------------------------- 116 pub fn cartesian_mesh_coords( 117 ceed: &Ceed, 118 dim: usize, 119 num_xyz: [usize; 3], 120 mesh_degree: usize, 121 mesh_size: usize, 122 ) -> libceed::Result<Vector> { 123 let p = mesh_degree + 1; 124 let mut num_d = [0; 3]; 125 let mut scalar_size = 1; 126 for d in 0..dim { 127 num_d[d] = num_xyz[d] * (p - 1) + 1; 128 scalar_size *= num_d[d]; 129 } 130 131 // Lobatto points 132 let lobatto_basis = 133 ceed.basis_tensor_H1_Lagrange(1, 1, 2, p, libceed::QuadMode::GaussLobatto)?; 134 let nodes_corners = ceed.vector_from_slice(&[0.0, 1.0])?; 135 let mut nodes_full = ceed.vector(p)?; 136 lobatto_basis.apply( 137 1, 138 libceed::TransposeMode::NoTranspose, 139 libceed::EvalMode::Interp, 140 &nodes_corners, 141 &mut nodes_full, 142 )?; 143 144 // Coordinates for mesh 145 let mut mesh_coords = ceed.vector(mesh_size)?; 146 mesh_coords.set_value(0.0)?; 147 { 148 let mut coords = mesh_coords.view_mut()?; 149 let nodes = nodes_full.view()?; 150 for gs_nodes in 0..scalar_size { 151 let mut r_nodes = gs_nodes; 152 for d in 0..dim { 153 let d_1d = r_nodes % num_d[d]; 154 coords[gs_nodes + scalar_size * d] = ((d_1d / (p - 1)) as libceed::Scalar 155 + nodes[d_1d % (p - 1)]) 156 / num_xyz[d] as libceed::Scalar; 157 r_nodes /= num_d[d]; 158 } 159 } 160 } 161 Ok(mesh_coords) 162 } 163 164 // ---------------------------------------------------------------------------- 165