diff --git a/README.md b/README.md index 289f422..6f75c98 100755 --- a/README.md +++ b/README.md @@ -44,3 +44,6 @@ Creating mesh and Caculating electrostatic potential distribution [fem](docs/sou

zhk_logo

+ +## Adaptive meshing +To use the latest features in gmsh, aka `Mesh.AutomaticMeshSizeField`, you have to uninstall your default gmsh and compile it from source using external `p4est` library. See instructions [gmsh_build](docs/gmsh_build.md) diff --git a/docs/gmsh_build.md b/docs/gmsh_build.md new file mode 100644 index 0000000..1464701 --- /dev/null +++ b/docs/gmsh_build.md @@ -0,0 +1,107 @@ +# How to build Gmsh with p4est and Hxt + +Regular `pip install gmsh` does NOT have p4est or Hxt. You have to build +Gmsh yourself from source to get them. + +## 1. Install tools + +```bash +brew install autoconf automake libtool open-mpi cmake fltk opencascade gmp +``` + +## 2. Build p4est + +```bash +cd ~ +git clone --recursive https://github.com/cburstedde/p4est.git +cd p4est +./bootstrap +mkdir -p build && cd build +../configure --enable-mpi --disable-shared CC=mpicc CXX=mpicxx --prefix=$HOME/local/p4est +make -j$(sysctl -n hw.ncpu) install +``` + +run `make check` to test your build - you should get `PASS` for all tests. +Your build should be installed in `$HOME/local/p4est`. Inside that directory you should find: +``` +$HOME/local/p4est/ +├── bin/ # any p4est utility executables, if built +├── include/ # p4est.h, p8est.h, sc.h, etc. +├── lib/ +│ ├── libp4est.a # static lib (since you used --disable-shared) +│ ├── libsc.a # sc is bundled as a dependency +│ └── pkgconfig/ +│ ├── p4est.pc +│ └── sc.pc +└── share/ # docs/misc, if any +``` + +## 3. Build Gmsh + +```bash +cd ~ +git clone https://gitlab.onelab.info/gmsh/gmsh.git +cd gmsh +mkdir build && cd build +cmake -DENABLE_P4EST=1 -DENABLE_HXT=1 -DENABLE_MPI=ON -DENABLE_OPENMP=1 \ + -DENABLE_BUILD_DYNAMIC=1 \ + -DCMAKE_PREFIX_PATH="$HOME/local/p4est;/opt/homebrew" \ + -DCMAKE_C_COMPILER=$(brew --prefix open-mpi)/bin/mpicc \ + -DCMAKE_CXX_COMPILER=$(brew --prefix open-mpi)/bin/mpicxx \ + -DJPEG_INCLUDE_DIR=/opt/homebrew/opt/jpeg/include \ + -DJPEG_LIBRARY=/opt/homebrew/opt/jpeg/lib/libjpeg.dylib \ + .. +make -j$(sysctl -n hw.ncpu) +make install +``` + +Your new gmsh binary is now at `/usr/local/bin/gmsh`. + +## 4. Check it worked + +```bash +/usr/local/bin/gmsh -info | grep -i p4est +/usr/local/bin/gmsh -info | grep -i hxt +``` + +Both should print builded app info. If empty, it didn't build in. + +## 5. One extra fix needed + +If you plan to use `-bgm`/`-size_field` (the adaptive meshing pipeline), +you need one code patch or it will crash. In +`src/mesh/automaticMeshSizeField.cpp`, find `sc_MPI_Init(&argc, &argv);` +near the top of `forestCreate()` and wrap it like this: + +```cpp +int already_init = 0; +MPI_Initialized(&already_init); +if (!already_init) { + mpiret = sc_MPI_Init(&argc, &argv); + SC_CHECK_MPI(mpiret); +} +``` + +Then rebuild: +```bash +cd ~/gmsh/build +make -j$(sysctl -n hw.ncpu) +make install +``` + +## 6. Make sure Python uses this build, not pip's + +```bash +python -c "import gmsh; print(gmsh.__file__)" +``` + +If it doesn't point to your new build, run: +```bash +pip uninstall gmsh +``` + +Add this to ~/.zshrc if you want it permanent: +```bash +echo 'export PYTHONPATH="/usr/local/lib:$PYTHONPATH"' >> ~/.zshrc +source ~/.zshrc +``` \ No newline at end of file diff --git a/zeroheliumkit/fem/gmsher.py b/zeroheliumkit/fem/gmsher.py index e7347ff..62904d9 100755 --- a/zeroheliumkit/fem/gmsher.py +++ b/zeroheliumkit/fem/gmsher.py @@ -152,10 +152,59 @@ def __post_init__(self): self.surface_currents = [] + +@dataclass +class AutomaticMeshSizeFieldSettings: + """ + Settings for the p4est/Hxt AutomaticMeshSizeField adaptive refinement pass. + + When enabled, a coarse reference mesh is generated first (controlled by + `ref_clscale`), then Gmsh's AutomaticMeshSizeField computes a p4est + octree size field from that reference mesh's curvature/feature data. + That field is then combined with any Box/Distance fields already + configured and used as the background field for the real mesh. + + Requires a Gmsh build with p4est enabled -- see the module-level NOTE + at the top of this file. If p4est isn't available, `gmsh.model.mesh.field.add` + will raise ("Unknown field type") when this is turned on. + + Args: + enabled (bool): turn the p4est/Hxt adaptive pass on. Defaults to False, + No behavior change for existing users of GMSHmaker. + ref_clscale (float): global element size scale used for the coarse + reference mesh. Smaller = finer reference mesh = more feature data for the size field to work + with, at the cost of a slower reference-mesh pass. Defaults to 0.2. + nPointsPerCircle (int): node density around curved features; higher + = finer resolution on curves. Defaults to 20. + nPointsPerGap (int): controls how finely narrow gaps between features + get resolved. Defaults to 25. + gradation (float): how fast element size is allowed to grow moving + away from a fine region. Close to 1.0 = gradual/smooth + transitions (more elements); higher = allows abrupt jumps to + coarser elements (fewer elements). Defaults to 1.1. + hBulk (float): target element size for open/background regions far + from any feature. -1 leaves it to Gmsh's automatic choice. + hMin (float): hard floor on element size. -1 leaves it automatic. + hMax (float): hard ceiling on element size. -1 leaves it automatic. + use_hxt_3d (bool): if True and meshing in 3D, sets + Mesh.Algorithm3D = 10 (Hxt) for the reference-mesh pass. + Defaults to True. + """ + enabled: bool = False + ref_clscale: float = 0.2 + nPointsPerCircle: int = 20 + nPointsPerGap: int = 25 + gradation: float = 1.1 + hBulk: float = -1 + hMin: float = -1 + hMax: float = -1 + use_hxt_3d: bool = True + @dataclass class MeshSettings: dim: int = 3 fields: dict = field(default_factory=dict) + automatic_mesh_size_field: AutomaticMeshSizeFieldSettings = field(default_factory=AutomaticMeshSizeFieldSettings) @dataclass @@ -810,6 +859,66 @@ def make_distance_threshold_field_mesh(self, distances: list[DistanceFieldMeshSe return field_ids + def make_automatic_mesh_size_field(self, dim: int) -> int: + """ + Runs the p4est-based AutomaticMeshSizeField pass and returns a field id + that can be combined with other fields (Box/Distance) via Min, same as + the rest of `setup_mesh_fields`. + + This does three things, all within the current Gmsh session (no CLI / + file round-trip needed): + 1. Generates a coarse reference mesh at `automatic_mesh_size_field.ref_clscale`, + since AutomaticMeshSizeField needs an existing mesh to compute + curvature/feature data from -- it cannot be evaluated on bare + geometry. + 2. Builds the p4est octree size field from that reference mesh. + 3. Clears the reference mesh so the field (not the reference mesh + itself) is what drives the real mesh generated later in + `create_mesh`. + + Requires a Gmsh build with p4est enabled -- see module-level NOTE. + + Args: + dim (int): mesh dimension (2 or 3) used for the reference pass. + + Returns: + int: id of the AutomaticMeshSizeField field. + """ + cfg = self.mesh.automatic_mesh_size_field + + prev_size_factor = gmsh.option.getNumber("Mesh.MeshSizeFactor") + prev_algo_3d = gmsh.option.getNumber("Mesh.Algorithm3D") if dim == 3 else None + + # 1. coarse reference mesh -- gives AutomaticMeshSizeField + # curvature/feature data to compute from + gmsh.option.setNumber("Mesh.MeshSizeFactor", cfg.ref_clscale) + if dim == 3 and cfg.use_hxt_3d: + gmsh.option.setNumber("Mesh.Algorithm3D", 10) # Hxt + gmsh.model.mesh.generate(dim) + + # 2. compute the p4est octree size field from the reference mesh + field_id = gmsh.model.mesh.field.add("AutomaticMeshSizeField") + gmsh.model.mesh.field.setNumber(field_id, "nPointsPerCircle", cfg.nPointsPerCircle) + gmsh.model.mesh.field.setNumber(field_id, "nPointsPerGap", cfg.nPointsPerGap) + gmsh.model.mesh.field.setNumber(field_id, "gradation", cfg.gradation) + if cfg.hBulk > 0: + gmsh.model.mesh.field.setNumber(field_id, "hBulk", cfg.hBulk) + if cfg.hMin > 0: + gmsh.model.mesh.field.setNumber(field_id, "hMin", cfg.hMin) + if cfg.hMax > 0: + gmsh.model.mesh.field.setNumber(field_id, "hMax", cfg.hMax) + + # 3. clear the reference mesh and restore mesh-size options so the + # following combined-field generate() call in create_mesh() + # is driven purely by the field(s), not leftover options + gmsh.model.mesh.clear() + gmsh.option.setNumber("Mesh.MeshSizeFactor", prev_size_factor) + if prev_algo_3d is not None: + gmsh.option.setNumber("Mesh.Algorithm3D", prev_algo_3d) + + return field_id + + def setup_mesh_fields(self): """ Build and set the background mesh field. @@ -830,12 +939,21 @@ def setup_mesh_fields(self): case _: print(f"Mesh field '{mesh_field_name}' is not recognized.") + if self.mesh.automatic_mesh_size_field and self.mesh.automatic_mesh_size_field.enabled: + field_ids.append(self.make_automatic_mesh_size_field(dim=self.mesh.dim)) + + if not field_ids: + return + minimum = gmsh.model.mesh.field.add("Min") gmsh.model.mesh.field.setNumbers(minimum, "FieldsList", field_ids) gmsh.model.mesh.field.setAsBackgroundMesh(minimum) gmsh.model.occ.synchronize() + + + def create_geo(self): path = self.save.with_suffix(".geo_unrolled") gmsh.write(str(path)) @@ -855,11 +973,19 @@ def create_mesh(self, dim=2): bar = alive_it([0], title='Gmsh generation ', length=3, spinner='elements', force_tty=True) try: for _ in bar: + + if dim == 3 and self.mesh.automatic_mesh_size_field.enabled and self.mesh.automatic_mesh_size_field.use_hxt_3d: + gmsh.option.setNumber("Mesh.Algorithm3D", 10) # Hxt + + else: + gmsh.option.setNumber("Mesh.Algorithm3D", 1) gmsh.model.mesh.generate(dim) + print("mesh is constructed") gmsh.model.mesh.setOrder(1) gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) gmsh.option.setNumber("Mesh.Binary", 0) + path = self.save.with_suffix(".msh") gmsh.write(str(path)) print("mesh saved")