How to properly get the mesh face on the top of a structure

[CalebDmArcher]

I am trying to only get the top surface of the entire structure by using mesh.facesTop but it returns the top surface of all layers. As shown below, it is a two-layer structure and the top surface area should be 0.001 but it returns 0.002.

import fipy as fp
from fipy import CellVariable, Gmsh3D, DiffusionTerm, Viewer, FaceVariable, ImplicitSourceTerm
from fipy.tools import numerix
import numpy as np
import matplotlib
matplotlib.use('TkAgg')

# Gmsh script as a string
gmsh_script = """
SetFactory("OpenCASCADE");

// Define layer dimensions
Lx = 0.05;
Ly = 0.02;
t_cu = 0.02;
t_fr4 = 0.02;

// Define mesh sizes for each layer
dx_cu = 0.001;
dx_fr4 = 0.001;

// Create Layer 1
Box(1) = {0, 0, 0, Lx, Ly, t_cu};

// Extrude Layer 2 from Layer 1
Extrude {0, 0, t_fr4} {
  Surface{6};
}

// Meshing properties for each layer
Field[1] = Box;
Field[1].VIn = dx_cu;
Field[1].VOut = dx_cu;
Field[1].XMin = 0; Field[1].XMax = Lx;
Field[1].YMin = 0; Field[1].YMax = Ly;
Field[1].ZMin = 0; Field[1].ZMax = t_cu;

Field[2] = Box;
Field[2].VIn = dx_fr4;
Field[2].VOut = dx_fr4;
Field[2].XMin = 0; Field[2].XMax = Lx;
Field[2].YMin = 0; Field[2].YMax = Ly;
Field[2].ZMin = t_cu; Field[2].ZMax = t_cu + t_fr4;

Field[3] = Min;
Field[3].FieldsList = {1,2};

Background Field = 3;

// Generate 3D mesh
Mesh 3;
"""

# Use Gmsh2D to create a mesh from the Gmsh script
mesh = Gmsh3D(gmsh_script)

top_faces = mesh.facesTop  # Identify the top faces
cell_areas = mesh._faceAreas[top_faces]  # Calculate the areas of the top faces
A = sum(cell_areas)  # Calculate the uniform heat flux q/A
print(f"Top surface area is {A} m^2")

So I modified the mesh part to 1 layer and it returns 0.001:

# Gmsh script as a string
gmsh_script = """
SetFactory("OpenCASCADE");

// Define layer dimensions
Lx = 0.05;
Ly = 0.02;
t_cu = 0.02;
t_fr4 = 0.02;

// Define mesh sizes for each layer
dx_cu = 0.001;
dx_fr4 = 0.001;

// Create Layer 1
Box(1) = {0, 0, 0, Lx, Ly, t_cu};

// Meshing properties for each layer
Field[1] = Box;
Field[1].VIn = dx_cu;
Field[1].VOut = dx_cu;
Field[1].XMin = 0; Field[1].XMax = Lx;
Field[1].YMin = 0; Field[1].YMax = Ly;
Field[1].ZMin = 0; Field[1].ZMax = t_cu;

Field[2] = Min;
Field[2].FieldsList = {1};

Background Field = 2;

// Generate 3D mesh
Mesh 2;

I also tried 3 layer which returns 0.003.
Since I need to mask only the top surface to do some calculations, may I know what is the proper way to only get the top surface?
(I also need to mask only the outer surfaces of the entire structure. I haven't tested yet, but I am afraid I might get the outer surfaces of all layers as well. )

[CalebDmArcher]

I just tried the code below to calculate the total surface area. Seems like it is correct.
Not sure why the top surface area is wrong.

exterior_faces = mesh.exteriorFaces
exterior_areas = mesh._faceAreas[exterior_faces]
total_exterior_area = sum(exterior_areas)
print(f"Total exterior surface area is {total_exterior_area} m^2")

The above code returns the result 0.0076 for the two-layer case.
(I thought it would return 0.0096 as it would count all outer surfaces including the overlapped face between the two layers, but it did not.)

[guyer]

I think you want mesh.facesFront or mesh.facesBack. Your boxes are stacked along the z direction. mesh.facesTop returns the faces with centers that maximize y.