Mesh Generating Tool Suggestion

Hi, are there any good mesh-generating tools suggested that are compatible with Fipy other than Gmsh?

I am using Gmsh currently but seems like Gmsh is not good at handling some cases such as a multiple-thin-layer structure: I used

mesh = Gmsh3D("cube.geo")
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)

to calculate the area of the top surface of the structure for the mesh file below. And get the value of A (3e-5) which is much smaller than if I calculate the area by using width times length (Lx Ly = 0.05 0.02 = 1e-4).

And there are also some ill-shape warnings which might imply that the mesh quality is not good. So even if I use (Lx * Ly) as the top surface area, the solution from solving the equation using Fipy is also not correct which I suspect is due to the poor mesh quality. (The reason that I can tell the result is not correct is because it does not match the result from a paper as well as my calculation by hand.)

A person from Gmsh forum said that Gmsh primarily generates isotropic instructed meshes so I feel like it is not very suitable for the case that the structure is thin and I need the mesh to have large scales on the x and y direction but a small scale on the z direction. I can set a smaller mesh size in general but it will make the program extremely slow.

SetFactory("OpenCASCADE");

// Define layer dimensions
Lx = 0.05;
Ly = 0.02;
t_cu = 5e-05;
t_fr4 = 0.0005;

// Define mesh sizes for each layer
dx_cu = 2.5e-04;
dx_fr4 = 2.5e-04;

// 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};
}

// Extrude Layer 3 from Layer 2
Extrude {0, 0, t_cu} {
  Surface{11};
}

// 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] = Box;
Field[3].VIn = dx_cu;
Field[3].VOut = dx_cu;
Field[3].XMin = 0; Field[3].XMax = Lx;
Field[3].YMin = 0; Field[3].YMax = Ly;
Field[3].ZMin = t_cu + t_fr4; Field[3].ZMax = t_cu + t_fr4 + t_cu;

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

Background Field = 4;

// Generate 3D mesh
Mesh 3;

Below is the code for solving the problem by using Fipy which might not be important for this discussion but just in case:

# Solver Setup
import os
# os.environ["FIPY_SOLVERS"] = "scipy"
os.environ["FIPY_SOLVERS"] = "petsc"
print("FIPY_SOLVERS is set to:", os.environ["FIPY_SOLVERS"])
os.environ["OMP_NUM_THREADS"] = "1"
print("OMP_NUM_THREADS is set to:", os.environ["OMP_NUM_THREADS"])

# Solver definitions
import fipy as fp
# scipy LU
# solver = fp.solvers.scipy.linearLUSolver.LinearLUSolver(tolerance=1e-10, iterations=10)

# scipy GMRES
# solver = fp.solvers.scipy.linearGMRESSolver.LinearGMRESSolver(iterations=1000, tolerance=1e-10)

# petsc LU without precon
# solver = fp.solvers.petsc.linearLUSolver.LinearLUSolver(tolerance=1e-10, iterations=10, precon='lu')

# petsc GMRES with precon, jacobi or bjacobi
solver = fp.solvers.petsc.linearGMRESSolver.LinearGMRESSolver(precon='jacobi', iterations=1000, tolerance=1e-10)

# petsc GMRES wihout precon
# solver = fp.solvers.petsc.linearGMRESSolver.LinearGMRESSolver(iterations=1000, tolerance=1e-10)

print("The solver package is being used is", fp.solvers.solver)

from fipy import CellVariable, Gmsh3D, DiffusionTerm, Viewer, FaceVariable, ImplicitSourceTerm
from fipy.tools import numerix
import numpy as np

import matplotlib
matplotlib.use('TkAgg')

import matplotlib.pyplot as plt
import time
from scipy.interpolate import griddata
# import gmsh

# Define layer dimensions
Lx = 5e-2
Ly = 2e-2
t_cu = 5e-5
t_fr4 = 5e-4

# Define mesh sizes for each layer and the small volume
dx_cu = 2.5e-5
dx_fr4 = 2.5e-5

# Define ambient temperature
T_ambient = 25.

# Define thermal conductivities
k_cu = 400.
k_fr4 = 0.5

# Define free air convection coefficient
h = 0.

# Define heat source
q = 50.

# # draw the mesh
# import gmsh
# # Initialize the Gmsh Python API
# gmsh.initialize()
# # gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2)
# # Merge the .geo file
# gmsh.merge('cube.geo')
# # Generate the 3D mesh
# gmsh.model.mesh.generate(3)
# # Launch the Gmsh GUI to visualize the mesh (optional)
# gmsh.fltk.run()
# # Finalize the Gmsh Python API
# gmsh.finalize()

mesh = Gmsh3D("cube.geo")

temp = CellVariable(name="solution variable", mesh=mesh, value=T_ambient)

# boundary conditions
# Bottom surface (fixed temperature)
bottom_faces = mesh.facesBottom
temp.constrain(T_ambient, bottom_faces)

# Side surfaces (adiabatic)

# Top surface (total heat flow)
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
# A = Lx * Ly
# q_A = q / A
# print(f"Top surface area is {A:.2f} m^2")
print(f"Top surface area is {A} m^2")
print(f"Heat power per unit area is {q_A:.2f} W/m^2")

# Set spatial variables for calculation and plot
xc, yc, zc = mesh.cellCenters
x, y, z = mesh.faceCenters

# Define the thermal conductivities for the equation
D = FaceVariable(mesh=mesh, value=k_cu)
mask_fr4 = (z > t_cu) & (z < (t_cu + t_fr4))
D.setValue(k_fr4, where=mask_fr4)

# Set the self-heating source
# No self-heating source

# Define the equation

# Tobin boundary condition
maskTopSurface = mesh.facesTop
Gamma = D * ~maskTopSurface
# Gamma = D
dPf = FaceVariable(mesh=mesh,
                   value=mesh._faceToCellDistanceRatio * mesh.cellDistanceVectors)

n = mesh.faceNormals
a = FaceVariable(mesh=mesh, value=h*n, rank=1)
b = D
g = FaceVariable(mesh=mesh, value= h*T_ambient + q_A, rank=0)

RobinCoeff = maskTopSurface * D * n / (dPf.dot(a) + b)

eq = (DiffusionTerm(coeff=Gamma) +
      # source +
      (RobinCoeff * g).divergence - ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence))

# start the timer
print('Start of the equation solving process...')
start_time = time.process_time()
start_time_e = time.perf_counter()

# Solve the equation
# # sweeping
# dx = dx_cu
# timeStepDuration = 0.9 * dx**2 / (2 * D) *10
# res = 1e+10
# previous_res = float('inf')
# tolerance = 1e-10
# sweeps = 0
#
# while abs(previous_res - res) > tolerance:
#     previous_res = res
#
#     # if solver is defined above (when using PETSc)
#     res = eq.sweep(var=temp, dt=timeStepDuration, solver=solver)
#
#     # if solver is not defined above
#     # res = eq.sweep(var=phi, dt=timeStepDuration)
#
#     sweeps += 1
#
#     print(f'Current residual is {res} of sweep number {sweeps}.')
#     print("Current sweep number is ", sweeps)

# noun sweeping
eq.solve(var=temp)

# end the timer
end_time_e = time.perf_counter()
end_time = time.process_time()
elapsed_time = end_time_e - start_time_e
cpu_time = end_time - start_time
print('//////////END of equation solving process//////////')
print(f"CPU time: {cpu_time:.4f}s; elapsed time: {elapsed_time:.4f}s.")
# print(f'converged at residual {res} of sweep number {sweeps-1}')

if __name__ == '__main__':
    # Extract the temperature values on the top surface
    top_surface_temps = temp.faceValue[top_faces]

    # Find and print the maximum temperature value on the top surface
    max_temp_top_surface = float(top_surface_temps.max())
    print(f"Maximum temperature on the top surface: {max_temp_top_surface:.2f} C")

    # Plot the temperature distribution on the top surface
    fig, ax = plt.subplots()
    sc = ax.scatter(x[top_faces], y[top_faces], c=numerix.array(top_surface_temps), cmap='hot')
    plt.colorbar(sc, ax=ax, label='Temperature (C)')
    ax.set_xlabel('x (m)')
    ax.set_ylabel('y (m)')
    ax.set_title('Temperature Distribution on Top Surface')
    plt.show()

usnistgov / fipy

Mesh Generating Tool Suggestion #1041