1using LibCEED, Printf, StaticArrays 2 3include("common.jl") 4 5function transform_mesh_coords!(dim, mesh_size, mesh_coords) 6 @witharray coords = mesh_coords begin 7 @inbounds @simd for i = 1:mesh_size 8 # map [0,1] to [0,1] varying the mesh density 9 coords[i] = 0.5 + 1.0/sqrt(3.0)*sin((2.0/3.0)*pi*(coords[i] - 0.5)) 10 end 11 end 12 exact_sa = (dim == 1 ? 2 : dim == 2 ? 4 : 6) 13end 14 15function run_ex2(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, gallery) 16 ncompx = dim 17 prob_size < 0 && (prob_size = 256*1024) 18 19 mesh_order = max(mesh_order, sol_order) 20 sol_order = mesh_order 21 22 ceed = Ceed(ceed_spec) 23 mesh_basis = 24 create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS) 25 sol_basis = 26 create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS) 27 28 nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size) 29 println("Mesh size: ", nxyz) 30 31 # Build CeedElemRestriction objects describing the mesh and solution discrete 32 # representations. 33 mesh_size, mesh_rstr, _ = build_cartesian_restriction( 34 ceed, 35 dim, 36 nxyz, 37 mesh_order, 38 ncompx, 39 num_qpts, 40 mode=RestrictionOnly, 41 ) 42 sol_size, _, qdata_rstr_i = build_cartesian_restriction( 43 ceed, 44 dim, 45 nxyz, 46 sol_order, 47 div(dim*(dim + 1), 2), 48 num_qpts, 49 mode=StridedOnly, 50 ) 51 sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction( 52 ceed, 53 dim, 54 nxyz, 55 sol_order, 56 1, 57 num_qpts, 58 mode=RestrictionAndStrided, 59 ) 60 println("Number of mesh nodes : ", div(mesh_size, dim)) 61 println("Number of solution nodes : ", sol_size) 62 63 # Create a CeedVector with the mesh coordinates. 64 mesh_coords = CeedVector(ceed, mesh_size) 65 set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords) 66 67 # Apply a transformation to the mesh. 68 exact_sa = transform_mesh_coords!(dim, mesh_size, mesh_coords) 69 70 # Create the Q-function that builds the diffusion operator (i.e. computes its 71 # quadrature data) and set its context data. 72 if !gallery 73 @interior_qf build_qfunc = ( 74 ceed, 75 dim=dim, 76 (J, :in, EVAL_GRAD, dim, dim), 77 (w, :in, EVAL_WEIGHT), 78 (qdata, :out, EVAL_NONE, dim*(dim + 1)÷2), 79 begin 80 Jinv = inv(J) 81 qdata .= setvoigt(w*det(J)*Jinv*Jinv') 82 end, 83 ) 84 else 85 build_qfunc = create_interior_qfunction(ceed, "Poisson$(dim)DBuild") 86 end 87 88 # Create the operator that builds the quadrature data for the diffusion 89 # operator. 90 build_oper = Operator( 91 ceed, 92 qf=build_qfunc, 93 fields=[ 94 (gallery ? :dx : :J, mesh_rstr, mesh_basis, CeedVectorActive()), 95 (gallery ? :weights : :w, ElemRestrictionNone(), mesh_basis, CeedVectorNone()), 96 (:qdata, qdata_rstr_i, BasisNone(), CeedVectorActive()), 97 ], 98 ) 99 100 # Compute the quadrature data for the diffusion operator. 101 elem_qpts = num_qpts^dim 102 num_elem = prod(nxyz) 103 qdata = CeedVector(ceed, num_elem*elem_qpts*div(dim*(dim + 1), 2)) 104 print("Computing the quadrature data for the diffusion operator ...") 105 flush(stdout) 106 apply!(build_oper, mesh_coords, qdata) 107 println(" done.") 108 109 # Create the Q-function that defines the action of the diffusion operator. 110 if !gallery 111 @interior_qf apply_qfunc = ( 112 ceed, 113 dim=dim, 114 (du, :in, EVAL_GRAD, dim), 115 (qdata, :in, EVAL_NONE, dim*(dim + 1)÷2), 116 (dv, :out, EVAL_GRAD, dim), 117 begin 118 dXdxdXdxT = getvoigt(qdata) 119 dv .= dXdxdXdxT*du 120 end, 121 ) 122 else 123 apply_qfunc = create_interior_qfunction(ceed, "Poisson$(dim)DApply") 124 end 125 126 # Create the diffusion operator. 127 oper = Operator( 128 ceed, 129 qf=apply_qfunc, 130 fields=[ 131 (:du, sol_rstr, sol_basis, CeedVectorActive()), 132 (:qdata, qdata_rstr_i, BasisNone(), qdata), 133 (:dv, sol_rstr, sol_basis, CeedVectorActive()), 134 ], 135 ) 136 137 # Compute the mesh surface area using the diff operator: 138 # sa = 1^T \cdot abs( K \cdot x). 139 print("Computing the mesh surface area using the formula: sa = 1^T.|K.x| ...") 140 flush(stdout) 141 142 # Create auxiliary solution-size vectors. 143 u = CeedVector(ceed, sol_size) 144 v = CeedVector(ceed, sol_size) 145 # Initialize 'u' with sum of coordinates, x+y+z. 146 u[] = 0.0 147 @witharray_read( 148 x_host = mesh_coords, 149 size = (mesh_size÷dim, dim), 150 @witharray(u_host = u, size = (sol_size, 1), sum!(u_host, x_host)) 151 ) 152 153 # Apply the diffusion operator: 'u' -> 'v'. 154 apply!(oper, u, v) 155 sa = witharray_read(x -> sum(abs, x), v, MEM_HOST) 156 157 println(" done.") 158 @printf("Exact mesh surface area : % .14g\n", exact_sa) 159 @printf("Computed mesh surface area : % .14g\n", sa) 160 @printf("Surface area error : % .14g\n", sa - exact_sa) 161end 162 163run_ex2( 164 ceed_spec="/cpu/self", 165 dim=3, 166 mesh_order=4, 167 sol_order=4, 168 num_qpts=6, 169 prob_size=-1, 170 gallery=false, 171) 172