ex2_Surface implementation in python

[arrowguy234]

Hello libCEED team

This PR is related to issue #1778. I'm happy to iterate on this further and open to any feedback or suggestions you have. Let me know if there are any style, structural, or interface adjustments you'd like to see!

Thanks a lot for your time

[jeremylt]

Ok, one general comment - I know its a bit tedious, but ex1 and ex2 are meant to explain to people how to use libCEED. In that spirit, can you add the comments from ex2-surface.c to this file? I know its a lot of comments, more than I'd usually recommend for a piece of software, but since this example is meant to teach the basic usage of libCEED, I think it is important.

[jeremylt]

CI is showing formatting issues. Please run make format-py in the root directory and commit the changes that autopep8 makes. Thanks

[arrowguy234]

@jeremylt I ran the format-py for formatting issues and also added comments similar to C code, I am little stuck on 2D and 3D area calculations In TransformMeshCoords function, is the error due to geometric mapping or approximation errors and should the exact surface area match exactly if the mesh transformation and operator setup are correct?

[jeremylt]

You won't be able to exactly replicate the surface area, but you should be able to reach the tolerances given here:

libCEED/examples/ceed/ex2-surface.c

Line 290 in 1095a7a

CeedScalar tol = (dim == 1 ? 10000. * CEED_EPSILON : dim == 2 ? 1E-1 : 1E-1);

[jeremylt]

The formatting test is showing these issues

--- a/examples/ceed/ex2_surface.py
+++ b/examples/ceed/ex2_surface.py
@@ -39,6 +39,8 @@ import argparse
 import libceed
 
 # Determine mesh size based on dimension, degree, and problem size
+
+
 def get_cartesian_mesh_size(dim, degree, prob_size):
     num_elem = prob_size // (degree ** dim)
     s = 0
@@ -56,6 +58,8 @@ def get_cartesian_mesh_size(dim, degree, prob_size):
     return num_xyz
 
 # Build Cartesian element restrictions
+
+
 def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts, for_qdata=False):
     p = degree + 1
     num_nodes = p ** dim
@@ -102,6 +106,8 @@ def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts,
     return restriction, size
 
 # Set mesh coordinates
+
+
 def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
     p = mesh_degree + 1
     nd = [num_xyz[d] * (p - 1) + 1 for d in range(dim)]
@@ -123,6 +129,8 @@ def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
     mesh_coords.set_array(coords, cmode=libceed.COPY_VALUES)
 
 # Transform mesh coordinates
+
+
 def transform_mesh_coords(dim, mesh_size, mesh_coords):
     with mesh_coords.array_write() as coords:
         for i in range(mesh_size):
@@ -130,6 +138,7 @@ def transform_mesh_coords(dim, mesh_size, mesh_coords):
     exact_surface = 2.0 if dim == 1 else 4.0 if dim == 2 else 6.0
     return exact_surface
 
+
 def main():
     # Process command line arguments
     parser = argparse.ArgumentParser()
@@ -156,7 +165,9 @@ def main():
     print(f"Approx. # unknowns     [-s] : {prob_size}")
 
     num_xyz = get_cartesian_mesh_size(dim, sol_degree, prob_size)
-    print(f"Mesh size: nx = {num_xyz[0]}" + (f", ny = {num_xyz[1]}" if dim > 1 else "") + (f", nz = {num_xyz[2]}" if dim > 2 else ""))
+    print(f"Mesh size: nx = {num_xyz[0]}" +
+          (f", ny = {num_xyz[1]}" if dim > 1 else "") +
+          (f", nz = {num_xyz[2]}" if dim > 2 else ""))
 
     mesh_restr, mesh_size = build_cartesian_restriction(ceed, dim, num_xyz, mesh_degree, dim, num_qpts)
     sol_restr, qdata_restr, sol_size, num_elem, elem_qpts = build_cartesian_restriction(
@@ -206,5 +217,6 @@ def main():
     print(f"Computed mesh surface area : {surface:.14g}")
     print(f"Surface area error         : {surface - exact_surface:.14g}")
 
+
 if __name__ == "__main__":
     main()

[arrowguy234]

@jeremylt


I think the code is giving acceptable mesh values according to the tolerances and its working for all 3 dimensions (1D, 2D, 3D)

[jeremylt]

Those are looking good. As long as you're matching the results from the C example, I think you can be confident this works as expected!

[jeremylt]

CI is reporting the following styling issues:

--- a/examples/ceed/ex2_surface.py
+++ b/examples/ceed/ex2_surface.py
@@ -38,6 +38,7 @@ import math
 import numpy as np
 import libceed
 
+
 def parse_arguments():
     parser = argparse.ArgumentParser(description="Compute curve length or surface area using libCEED")
     parser.add_argument("-c", "--ceed", default="/cpu/self",
@@ -58,6 +59,7 @@ def parse_arguments():
                         help="Use gallery QFunction")
     return parser.parse_args()
 
+
 def get_cartesian_mesh_size(dim, degree, prob_size):
     num_elem = prob_size // (degree ** dim)
     s = 0
@@ -75,6 +77,7 @@ def get_cartesian_mesh_size(dim, degree, prob_size):
         num_xyz.append(1 << sd)
     return num_xyz
 
+
 def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts):
     p = degree + 1
     num_nodes = p ** dim
@@ -108,7 +111,7 @@ def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts):
                 r_node //= p
             elem_nodes[e * num_nodes + n] = g_node
 
-    elem_restriction = ceed.ElemRestriction(num_elem, num_nodes, num_comp, 
+    elem_restriction = ceed.ElemRestriction(num_elem, num_nodes, num_comp,
                                             scalar_size, size, elem_nodes)
 
     qd_comp = dim * (dim + 1) // 2
@@ -118,6 +121,7 @@ def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts):
                                                      strides)
     return elem_restriction, size, q_data_restriction, num_elem, elem_qpts
 
+
 def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
     p = mesh_degree + 1
     nd = []
@@ -141,15 +145,18 @@ def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
 
     mesh_coords.set_array(coords, cmode=libceed.COPY_VALUES)
 
+
 def transform_mesh_coords(dim, mesh_size, mesh_coords):
     exact_measure = {1: 2.0, 2: 4.0, 3: 6.0}[dim]
     with mesh_coords.array() as coords:
         for d in range(dim):
             offset = d * (mesh_size // dim)
             for i in range(mesh_size // dim):
-                coords[offset + i] = 0.5 + (1.0/np.sqrt(3.0)) * np.sin((2.0/3.0) * np.pi * (coords[offset + i] - 0.5))
+                coords[offset + i] = 0.5 + (1.0 / np.sqrt(3.0)) * np.sin((2.0 / 3.0)
+                                                                         * np.pi * (coords[offset + i] - 0.5))
     return exact_measure
 
+
 def main():
     args = parse_arguments()
     ceed = libceed.Ceed(args.ceed)
@@ -169,10 +176,10 @@ def main():
     print(f"  Approx problem size: {args.problem_size}")
 
     num_xyz = get_cartesian_mesh_size(dim, sol_degree, args.problem_size)
-    print(f"Mesh size: {', '.join(f'n{chr(120+d)} = {num_xyz[d]}' for d in range(dim))}")
+    print(f"Mesh size: {', '.join(f'n{chr(120 + d)} = {num_xyz[d]}' for d in range(dim))}")
 
-    mesh_basis = ceed.BasisTensorH1Lagrange(dim, num_comp_x, mesh_degree+1, num_qpts, libceed.GAUSS)
-    sol_basis = ceed.BasisTensorH1Lagrange(dim, 1, sol_degree+1, num_qpts, libceed.GAUSS)
+    mesh_basis = ceed.BasisTensorH1Lagrange(dim, num_comp_x, mesh_degree + 1, num_qpts, libceed.GAUSS)
+    sol_basis = ceed.BasisTensorH1Lagrange(dim, 1, sol_degree + 1, num_qpts, libceed.GAUSS)
 
     mesh_restr, mesh_size, _, _, _ = build_cartesian_restriction(
         ceed, dim, num_xyz, mesh_degree, num_comp_x, num_qpts)
@@ -225,5 +232,6 @@ def main():
 
     return 0
 
+
 if __name__ == "__main__":
     sys.exit(main())

[arrowguy234]

@jeremylt The formatting changes have been completed. Please review the updated file and let me know if any further modifications are needed.

[jeremylt]

I'm still seeing these formatting issues:

diff --git a/examples/ceed/ex2_surface.py b/examples/ceed/ex2_surface.py
index bd4ade8..dbf81dc 100644
--- a/examples/ceed/ex2_surface.py
+++ b/examples/ceed/ex2_surface.py
@@ -45,6 +45,7 @@ import math
 import numpy as np
 import libceed
 
+
 def parse_arguments():
     parser = argparse.ArgumentParser(description="Compute curve length or surface area using libCEED")
     parser.add_argument("-c", "--ceed", default="/cpu/self",
@@ -65,6 +66,7 @@ def parse_arguments():
                         help="Use gallery QFunction")
     return parser.parse_args()
 
+
 def get_cartesian_mesh_size(dim, degree, prob_size):
     num_elem = prob_size // (degree ** dim)
     s = 0
@@ -81,6 +83,7 @@ def get_cartesian_mesh_size(dim, degree, prob_size):
         num_xyz.append(1 << sd)
     return num_xyz
 
+
 def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts):
     p = degree + 1
     num_nodes = p ** dim
@@ -117,6 +120,7 @@ def build_cartesian_restriction(ceed, dim, num_xyz, degree, num_comp, num_qpts):
         num_elem, elem_qpts, qd_comp, num_elem * elem_qpts * qd_comp, strides)
     return elem_restriction, size, q_data_restriction, num_elem, elem_qpts
 
+
 def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
     p = mesh_degree + 1
     nd = []
@@ -137,16 +141,18 @@ def set_cartesian_mesh_coords(ceed, dim, num_xyz, mesh_degree, mesh_coords):
             r_node //= nd[d]
     mesh_coords.set_array(coords, cmode=libceed.COPY_VALUES)
 
+
 def transform_mesh_coords(dim, mesh_size, mesh_coords):
     exact_measure = {1: 2.0, 2: 4.0, 3: 6.0}[dim]
     with mesh_coords.array() as coords:
         for d in range(dim):
             offset = d * (mesh_size // dim)
             for i in range(mesh_size // dim):
-                coords[offset + i] = 0.5 + (1.0/np.sqrt(3.0)) * np.sin(
-                    (2.0/3.0) * np.pi * (coords[offset + i] - 0.5))
+                coords[offset + i] = 0.5 + (1.0 / np.sqrt(3.0)) * np.sin(
+                    (2.0 / 3.0) * np.pi * (coords[offset + i] - 0.5))
     return exact_measure
 
+
 def main():
     args = parse_arguments()
     ceed = libceed.Ceed(args.ceed)
@@ -162,11 +168,11 @@ def main():
     print(f"  Quadrature points: {num_qpts}")
     print(f"  Approx problem size: {args.problem_size}")
     num_xyz = get_cartesian_mesh_size(dim, sol_degree, args.problem_size)
-    print(f"Mesh size: {', '.join(f'n{chr(120+d)} = {num_xyz[d]}' for d in range(dim))}")
+    print(f"Mesh size: {', '.join(f'n{chr(120 + d)} = {num_xyz[d]}' for d in range(dim))}")
     mesh_basis = ceed.BasisTensorH1Lagrange(
-        dim, num_comp_x, mesh_degree+1, num_qpts, libceed.GAUSS)
+        dim, num_comp_x, mesh_degree + 1, num_qpts, libceed.GAUSS)
     sol_basis = ceed.BasisTensorH1Lagrange(
-        dim, 1, sol_degree+1, num_qpts, libceed.GAUSS)
+        dim, 1, sol_degree + 1, num_qpts, libceed.GAUSS)
     mesh_restr, mesh_size, _, _, _ = build_cartesian_restriction(
         ceed, dim, num_xyz, mesh_degree, num_comp_x, num_qpts)
     sol_restr, sol_size, q_data_restr, num_elem, elem_qpts = build_cartesian_restriction(
@@ -208,5 +214,6 @@ def main():
     print(f"Error:                  {measure - exact_measure:.3e}")
     return 0
 
+
 if __name__ == "__main__":
     sys.exit(main()

[jeremylt]

Suggested change

	def main():
	args = parse_arguments()
	def example_2(args):

See the discussion here: #1811

You'll want a main below this fn with the contents

def main():
    """Main function for surface area computation example"""
    args = parse_arguments()
    return example_2(args)

[arrowguy234]

I’ve reviewed the discussion and updated my code to follow the structure used in the volume example, including adjusting object usage and function layout. However, I’m still having trouble with formatting.
I ran make format-py, and while it shows formatting changes in the diff, those changes aren't being committed to the branch. Could you help me ensure that the formatted code is properly staged and committed?

[jeremylt]

The comments and update to main look good - thanks

For the formatting, we use autopep8. You could try to use it directly?

$ autopep8 --in-place --aggressive --max-line-length 120

then

$ git add .

then

$ git commit -m "style fixes"

and push

[arrowguy234]

@jeremylt I used the autopep8 command directly, and I believe all formatting changes have now been applied. Please have a look and let me know if anything else needs to be updated.

[jeremylt]

This looks pretty close to me - I added a handful of comments for small changes

[arrowguy234]

@jeremylt I believe all requested changes have been addressed, and the file has been moved under examples/python. Please let me know if there are any further updates needed.

[jeremylt]

This looks good to me. Thanks for the hard work!

PR accepted (CC: @valeriabarra)

Let me know when you're ready for me to merge this to the staging branch

[arrowguy234]

@jeremylt Thanks so much! I'm glad the changes look good.Really appreciate your detailed feedback and support throughout the process!

[valeriabarra]

Nice job @arrowguy234 ! I am glad this contribution can be merged in!

[arrowguy234]

@valeriabarra Thanks a lot for giving me this opportunity.