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