xref: /libCEED/examples/python/ex2_surface.py (revision d4cc18453651bd0f94c1a2e078b2646a92dafdcc) 
1#!/usr/bin/env python3
2# Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
3# All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
4#
5# SPDX-License-Identifier: BSD-2-Clause
6#
7# This file is part of CEED:  http://github.com/ceed
8#
9# libCEED example using diffusion operator to compute surface area
10#
11# Sample runs:
12#
13#     python ex2_surface.py
14#     python ex2_surface.py -c /cpu/self
15#     python ex2_surface.py -c /gpu/cuda
16
17import sys
18import os
19import numpy as np
20import libceed
21import ex_common as common
22
23
24def main():
25    """Main driver for surface area example"""
26    args = common.parse_arguments()
27    return example_2(args)
28
29
30def example_2(options):
31    """Compute surface area using diffusion operator
32
33    Args:
34        args: Parsed command line arguments
35
36    Returns:
37        int: 0 on success, error code on failure
38    """
39    # Process arguments
40    args = options
41    dim = args.dim
42    mesh_degree = max(args.mesh_degree, args.solution_degree)
43    sol_degree = args.solution_degree
44    num_qpts = args.quadrature_points
45    problem_size = args.problem_size if args.problem_size > 0 else (500 * dim * dim if args.test else 256 * 1024)
46    ncomp_x = dim  # Number of coordinate components
47
48    # Print configuration
49    if not args.quiet:
50        print("Selected options: [command line option] : <current value>")
51        print(f"    Ceed specification [-c] : {args.ceed}")
52        print(f"    Mesh dimension     [-d] : {dim}")
53        print(f"    Mesh degree        [-m] : {mesh_degree}")
54        print(f"    Solution degree    [-p] : {sol_degree}")
55        print(f"    Num. 1D quadr. pts [-q] : {num_qpts}")
56        print(f"    Approx. # unknowns [-s] : {problem_size}")
57        print(f"    QFunction source   [-g] : {'gallery' if args.gallery else 'user'}")
58
59    # Initialize CEED
60    ceed = libceed.Ceed(args.ceed)
61
62    # Create bases
63    # Tensor-product Lagrange basis for mesh coordinates
64    mesh_basis = ceed.BasisTensorH1Lagrange(
65        dim, ncomp_x, mesh_degree + 1, num_qpts, libceed.GAUSS)
66
67    # Tensor-product Lagrange basis for solution
68    solution_basis = ceed.BasisTensorH1Lagrange(
69        dim, 1, sol_degree + 1, num_qpts, libceed.GAUSS)
70
71    # Create mesh
72    # Determine mesh size
73    num_xyz = common.get_cartesian_mesh_size(dim, sol_degree, problem_size)
74    if not args.quiet:
75        print("\nMesh size                   : nx = %d" % num_xyz[0], end="")
76        if dim > 1:
77            print(", ny = %d" % num_xyz[1], end="")
78        if dim > 2:
79            print(", nz = %d" % num_xyz[2], end="")
80        print()
81
82    # Create element restrictions
83    num_q_comp = dim * (dim + 1) // 2
84    mesh_restriction, mesh_size, _, _, _ = common.build_cartesian_restriction(
85        ceed, dim, num_xyz, mesh_degree, ncomp_x, num_q_comp, num_qpts, create_qdata=False)
86    solution_restriction, sol_size, q_data_restriction, num_elem, elem_qpts = common.build_cartesian_restriction(
87        ceed, dim, num_xyz, sol_degree, 1, num_q_comp, num_qpts, create_qdata=True)
88
89    if not args.quiet:
90        print("Number of mesh nodes        : %d" % (mesh_size // dim))
91        print("Number of solution nodes    : %d" % sol_size)
92
93    # Create and transform mesh coordinates
94    mesh_coords = ceed.Vector(mesh_size)
95    common.set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords)
96    _, exact_surface_area = common.transform_mesh_coords(dim, mesh_size, mesh_coords, use_sin=False)
97
98    # Create the QFunction that builds the diffusion operator (i.e. computes
99    # its quadrature data) and set its context data
100    qf_build = None
101    if args.gallery:
102        qf_build = ceed.QFunctionByName(f"Poisson{dim}DBuild")
103    else:
104        build_ctx = ceed.QFunctionContext()
105        ctx_data = np.array([dim, dim], dtype=np.int32)
106        build_ctx.set_data(ctx_data)
107
108        qfs_so = common.load_qfs_so()
109        file_dir = os.path.dirname(os.path.abspath(__file__))
110
111        qf_build = ceed.QFunction(1, qfs_so.build_diff,
112                                  os.path.join(file_dir, "ex2-surface.h:build_diff"))
113        qf_build.add_input("dx", dim * dim, libceed.EVAL_GRAD)
114        qf_build.add_input("weights", 1, libceed.EVAL_WEIGHT)
115        qf_build.add_output("qdata", num_q_comp, libceed.EVAL_NONE)
116        qf_build.set_context(build_ctx)
117
118    # Operator for building quadrature data
119    op_build = ceed.Operator(qf_build)
120    op_build.set_field("dx", mesh_restriction, mesh_basis, libceed.VECTOR_ACTIVE)
121    op_build.set_field("weights", libceed.ELEMRESTRICTION_NONE, mesh_basis, libceed.VECTOR_NONE)
122    op_build.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, libceed.VECTOR_ACTIVE)
123
124    # Compute quadrature data
125    q_data = ceed.Vector(num_elem * elem_qpts * num_q_comp)
126    op_build.apply(mesh_coords, q_data)
127
128    # Create the QFunction that defines the action of the diffusion operator
129    qf_diff = None
130    if args.gallery:
131        qf_diff = ceed.QFunctionByName(f"Poisson{dim}DApply")
132    else:
133        build_ctx = ceed.QFunctionContext()
134        ctx_data = np.array([dim, dim], dtype=np.int32)
135        build_ctx.set_data(ctx_data)
136
137        qfs_so = common.load_qfs_so()
138        file_dir = os.path.dirname(os.path.abspath(__file__))
139
140        qf_diff = ceed.QFunction(1, qfs_so.apply_diff,
141                                 os.path.join(file_dir, "ex2-surface.h:apply_diff"))
142        qf_diff.add_input("du", dim, libceed.EVAL_GRAD)
143        qf_diff.add_input("qdata", num_q_comp, libceed.EVAL_NONE)
144        qf_diff.add_output("dv", dim, libceed.EVAL_GRAD)
145        qf_diff.set_context(build_ctx)
146
147    # Diffusion operator
148    op_diff = ceed.Operator(qf_diff)
149    op_diff.set_field("du", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE)
150    op_diff.set_field("qdata", q_data_restriction, libceed.BASIS_NONE, q_data)
151    op_diff.set_field("dv", solution_restriction, solution_basis, libceed.VECTOR_ACTIVE)
152
153    # Create vectors
154    u = ceed.Vector(sol_size)  # Input vector
155    v = ceed.Vector(sol_size)  # Output vector
156
157    # Initialize u with sum of coordinates (x + y + z)
158    with mesh_coords.array_read() as x_array, u.array_write() as u_array:
159        for i in range(sol_size):
160            u_array[i] = sum(x_array[i + j * (sol_size)] for j in range(dim))
161
162    # Apply operator: v = K * u
163    op_diff.apply(u, v)
164
165    # Compute surface area by summing absolute values of v
166    surface_area = 0.0
167    with v.array_read() as v_array:
168        surface_area = np.sum(abs(v_array))
169
170    if not args.test:
171        print()
172        print(f"Exact mesh surface area    : {exact_surface_area:.14g}")
173        print(f"Computed mesh surface area : {surface_area:.14g}")
174        print(f"Surface area error         : {surface_area - exact_surface_area:.14g}")
175    else:
176        # Test mode - check if error is within tolerance
177        tol = 10000 * libceed.EPSILON if dim == 1 else 1e-1
178        if abs(surface_area - exact_surface_area) > tol:
179            print(f"Surface area error : {surface_area - exact_surface_area:.14g}")
180            sys.exit(1)
181
182    return 0
183
184
185if __name__ == "__main__":
186    sys.exit(main())
187