diff --git a/.github/workflows/.readthedocs.yaml b/.readthedocs.yaml similarity index 100% rename from .github/workflows/.readthedocs.yaml rename to .readthedocs.yaml diff --git a/CMakePresets.json b/CMakePresets.json index e30fe15..d3928e6 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -46,6 +46,17 @@ "MFEM_ENABLE_TESTING": false, "MFEM_USE_OPENMP": false } + }, + { + "name": "gnu-dbg", + "displayName": "GNU g++ compiler", + "inherits": "default", + "cacheVariables": { + "TULIP_USE_MFEM_AS_SUBDIRECTORY": true, + "MFEM_ENABLE_TESTING": false, + "MFEM_USE_OPENMP": false, + "CMAKE_BUILD_TYPE": "RelWithDebInfo" + } } ], "buildPresets": [ diff --git a/docs/tulip_data_format.md b/docs/tulip_data_format.md index 5c2f1d7..947d058 100644 --- a/docs/tulip_data_format.md +++ b/docs/tulip_data_format.md @@ -7,7 +7,7 @@ - [`shield`](#shield) - [`dielectric`](#dielectric) - [`open`](#open) - - [``](#model) + - [``](#layers) - [.tulip.adapted.json file format](#tulipadaptedjson-file-format) - [Example](#example) - [`materials` array](#materials-array) @@ -21,12 +21,14 @@ Tulip uses three types of file formats: - `CASE_NAME.tulip.out.json` which is the solver output containing the $L$ and $C$ PUL matrices for shielded domains and the multipolar expansion coefficients for an open domain. # .tulip.input.json file format -Tulip receives a JSON object as an input with the entries described below. Square brackets indicate that the entry is optional and a default value will be assumed, angle brackets indicate that the entry is mandatory. +Tulip receives a JSON object as an input with the entries described below. Square brackets indicate that the entry is optional and a default value will be assumed, angle brackets indicate that the entry is mandatory. Unless specified otherwise all units are assumed to be in SI-MKS. Filename should be in the format `CASE_NAME.tulip.input.json`. +At minimum, the input JSON must include top-level `materials` and `layers` arrays. + ## `[adapterOptions]` It can contain the following entries, as explained in [AdapterOption.h](../src/adapter/AdapterOptions.h) with their corresponding default values. An example is shown below. ```json @@ -49,7 +51,7 @@ Driver manages the solver`and generates outputs. Default options can be checked ``` ## `` -These materials are associated with `model` `layers` to define regions with different material properties. +These materials are associated with top-level `layers` to define regions with different material properties. They are defined by an array of JSON objects with: - `[name]` a string with a human readable name. - `` an integer identifier with a unique number. @@ -84,12 +86,11 @@ A dielectric is defined with a `[relativePermittivity]` which defaults to `1.0`. An `open` material serves to specify the computational boundary of the problem. It must intersect every other material layer. If no open boundary is specified for an open problem, one is computed automatically, together with _inner_ and _outer_ regions used to extract the unshielded multiwire coefficients. -## `` -This object can contain the following entries: -+ `` which is an array which associates the layers present in the `.step` file with the different `materials`. Each layer is specified by: - - `` which must match exactly the name of the corresponding layer within the `.step` file. It must be unique. - - `` which is an integer non-negative unique identifier which will be used to order the results for the calculated PUL matrices. - - `` which must match an `id` from a material in the list of `materials` +## `` +This top-level array associates the layers present in the `.step` file with the different `materials`. Each layer is specified by: +- `` which must match exactly the name of the corresponding layer within the `.step` file. It must be unique. +- `` which is an integer non-negative unique identifier which will be used to order the results for the calculated PUL matrices. +- `` which must match an `id` from a material in the list of `materials` # .tulip.adapted.json file format @@ -177,4 +178,4 @@ This object can contain the following entries: For unshielded-domains stores the parameters needed to reconstruct the field using a multipolar expansion. - It also stores `materialAssociation` information which serves to reconstruct the \ No newline at end of file + It also stores `materialAssociation` information which serves to reconstruct the diff --git a/scripts/compute_incell_parameters.py b/scripts/compute_incell_parameters.py new file mode 100644 index 0000000..13a0ac3 --- /dev/null +++ b/scripts/compute_incell_parameters.py @@ -0,0 +1,329 @@ +""" +compute_incell_parameters.py +============================ +Postprocessing script for a Tulip solved case. + +Reads a .tulip.out.json file and, for each material with a multipolar +expansion, computes the in-cell capacitance C[i,j] (F/m) and inductance +L[i,j] (H/m) for a fixed reference conductor i and all available conductors j +in that material association. + +Formulas (from InCellPotentials::getInCellCapacitanceUsingInnerRegion and +getInCellInductanceUsingInnerRegion in Results.cpp): + + C[i,j] = Q_j / (V_i|Vj=1 - _inner) * epsilon0 + + L[i,j] = (A_i|Aj=1 - _inner) / I_j * mu0 + +where: + Q_j = electric[j].ab[0][0] (monopole charge coefficient) + = electric[j].innerRegionAveragePotential + V_i = electric[j].conductorPotentials[i] + + I_j = magnetic[j].ab[0][0] (monopole current coefficient) + = magnetic[j].innerRegionAveragePotential + A_i = magnetic[j].conductorPotentials[i] + +Usage: + python scripts/compute_incell_parameters.py + python scripts/compute_incell_parameters.py --run-tests +""" + +import argparse +import contextlib +import io +import json +import math +import os +import sys +import unittest + +# Physical constants (match tulip/src/driver/constants.h) +EPSILON0_SI = 8.8541878176e-12 # F/m +MU0_SI = 4.0e-7 * math.pi # H/m + + +def load_json(json_path: str) -> dict: + """Load a Tulip output JSON file from disk.""" + abs_path = os.path.abspath(json_path) + with open(abs_path, encoding="utf-8") as fh: + return json.load(fh) + + +def capacitance(electric_solutions: list, i: int, j: int) -> float: + """Compute C[i,j] in F/m.""" + sol_j = electric_solutions[j] + q_j = sol_j["ab"][0][0] + avg_v_j = sol_j["innerRegionAveragePotential"] + v_i = sol_j["conductorPotentials"][i] + return q_j / (v_i - avg_v_j) * EPSILON0_SI + + +def inductance(magnetic_solutions: list, i: int, j: int) -> float: + """Compute L[i,j] in H/m.""" + sol_j = magnetic_solutions[j] + i_j = sol_j["ab"][0][0] + avg_a_j = sol_j["innerRegionAveragePotential"] + a_i = sol_j["conductorPotentials"][i] + return (a_i - avg_a_j) / i_j * MU0_SI + + +def local_index_from_element_id(element_ids: list, element_id: int): + """Map a physical conductor id to the local multipolar solution index.""" + try: + return element_ids.index(element_id) + except ValueError: + return None + + +def element_ids_for_material(material_id: int, mat_assoc: list, fallback_size: int) -> list: + """Return material-associated conductor ids, or fallback to local indices.""" + assoc = next((a for a in mat_assoc if a.get("materialId") == material_id), None) + if assoc is not None and "elementIds" in assoc: + return assoc["elementIds"] + return list(range(fallback_size)) + + +def compute_material_rows(material: dict, mat_assoc: list, i_ref: int): + """Compute (j, C[i_ref,j], L[i_ref,j]) rows for one material.""" + mat_id = material.get("id") + mat_type = material.get("type", "unknown") + mp = material.get("multipolarExpansion", {}) + e_sols = mp.get("electric", []) + m_sols = mp.get("magnetic", []) + + if not e_sols or not m_sols: + return None, "missing electric or magnetic solutions" + + element_ids = element_ids_for_material(mat_id, mat_assoc, len(e_sols)) + i_local = local_index_from_element_id(element_ids, i_ref) + if i_local is None: + return None, f"reference conductor id {i_ref} not in {element_ids}" + + rows = [] + for j in element_ids: + j_local = local_index_from_element_id(element_ids, j) + if j_local is None: + continue + c_val = capacitance(e_sols, i_local, j_local) + l_val = inductance(m_sols, i_local, j_local) + rows.append((j, c_val, l_val)) + + return { + "mat_id": mat_id, + "mat_type": mat_type, + "element_ids": element_ids, + "rows": rows, + }, None + + +def print_material_report(report: dict, i_ref: int): + """Print computed rows for one material.""" + print(f"\n{'=' * 55}") + print(f" Material id={report['mat_id']} type={report['mat_type']}") + print(f" Total conductors: {len(report['element_ids'])}") + print(f"{'=' * 55}") + print(f" {'j':>4} {f'C[{i_ref},j] (F/m)':>18} {f'L[{i_ref},j] (H/m)':>18}") + print(f" {'-' * 4} {'-' * 18} {'-' * 18}") + + for j, c_val, l_val in report["rows"]: + print(f" {j:>4} {c_val:>+18.6e} {l_val:>+18.6e}") + + +def parse_args(argv: list): + """Parse command-line arguments.""" + parser = argparse.ArgumentParser( + description="Compute in-cell C/L parameters from a .tulip.out.json file." + ) + parser.add_argument( + "--run-tests", + action="store_true", + help="Run the unit tests embedded in this file.", + ) + parser.add_argument( + "json_path", + nargs="?", + help="Path to a .tulip.out.json file.", + ) + parser.add_argument( + "i_ref", + nargs="?", + type=int, + help="Reference conductor id i.", + ) + + args = parser.parse_args(argv) + if args.run_tests: + return args + + if args.json_path is None or args.i_ref is None: + parser.error("json_path and i_ref are required unless --run-tests is used") + return args + + +def run_self_tests() -> int: + """Execute same-file unit tests.""" + test_program = unittest.main(module=__name__, argv=[sys.argv[0]], exit=False) + return 0 if test_program.result.wasSuccessful() else 1 + + +def run(argv=None) -> int: + """Entrypoint used by both CLI and tests.""" + args = parse_args(sys.argv[1:] if argv is None else argv) + if args.run_tests: + return run_self_tests() + + json_path = os.path.abspath(args.json_path) + i_ref = args.i_ref + + print(f"Output JSON path : {json_path}") + print(f"Reference conductor i : {i_ref}") + + data = load_json(json_path) + materials = data.get("materials", []) + mat_assoc = data.get("materialAssociations", []) + + if not materials: + print("No materials found in output JSON.") + return 0 + + for material in materials: + report, warn = compute_material_rows(material, mat_assoc, i_ref) + if warn is not None: + print(f"[warn] Material {material.get('id')}: {warn}, skipping.") + continue + print_material_report(report, i_ref) + print() + + return 0 + + +class TestComputeInCellParameters(unittest.TestCase): + """Unit tests for compute_incell_parameters.py.""" + + def test_local_index_lookup(self): + element_ids = [16, 30, 42] + self.assertEqual(local_index_from_element_id(element_ids, 30), 1) + self.assertIsNone(local_index_from_element_id(element_ids, 99)) + + def test_compute_material_rows_uses_all_conductors(self): + material = { + "id": 7, + "type": "dielectric", + "multipolarExpansion": { + "electric": [ + { + "ab": [[1.0]], + "innerRegionAveragePotential": 0.0, + "conductorPotentials": [2.0, 3.0], + }, + { + "ab": [[2.0]], + "innerRegionAveragePotential": 1.0, + "conductorPotentials": [3.0, 4.0], + }, + ], + "magnetic": [ + { + "ab": [[1.0]], + "innerRegionAveragePotential": 0.0, + "conductorPotentials": [2.0, 3.0], + }, + { + "ab": [[2.0]], + "innerRegionAveragePotential": 1.0, + "conductorPotentials": [3.0, 4.0], + }, + ], + }, + } + mat_assoc = [{"materialId": 7, "elementIds": [16, 30]}] + + report, warn = compute_material_rows(material, mat_assoc, 16) + self.assertIsNone(warn) + self.assertIsNotNone(report) + self.assertEqual([row[0] for row in report["rows"]], [16, 30]) + + def test_reference_case_values_from_embedded_fixture(self): + material = { + "id": 1, + "type": "unshieldedMultiwire", + "multipolarExpansion": { + "electric": [ + { + "ab": [ + [0.8419073593680602, 0.0], + [4.7094817324223296e-07, -7.437578541169077e-07], + [-2.254972536944232e-06, -1.3845721768518446e-06], + [-1.7502712357816652e-09, 1.6065198400637674e-09], + ], + "conductorPotentials": [1.0, 0.7525321728251249], + "expansionCenter": [0.0017138274803646352, -0.0027066091802577024], + "innerRegionAveragePotential": 0.6100537777863153, + }, + { + "ab": [ + [1.0308546965043754, 0.0], + [4.236110742075576e-07, -8.85263389452124e-07], + [-1.301845581534472e-07, -8.009832224471683e-08], + [3.712118282384746e-10, -3.4152957909780026e-10], + ], + "conductorPotentials": [0.9213984270827741, 1.0], + "expansionCenter": [0.0023197106176927533, -0.004847736541849076], + "innerRegionAveragePotential": 0.7447102973543827, + }, + ], + "magnetic": [ + { + "ab": [ + [0.8419073593680602, 0.0], + [4.7094817324223296e-07, -7.437578541169077e-07], + [-2.254972536944232e-06, -1.3845721768518446e-06], + [-1.7502712357816652e-09, 1.6065198400637674e-09], + ], + "conductorPotentials": [1.0, 0.7525321728251249], + "expansionCenter": [0.0017138274803646352, -0.0027066091802577024], + "innerRegionAveragePotential": 0.6100537777863153, + }, + { + "ab": [ + [1.0308546965043754, 0.0], + [4.236110742075576e-07, -8.85263389452124e-07], + [-1.301845581534472e-07, -8.009832224471683e-08], + [3.712118282384746e-10, -3.4152957909780026e-10], + ], + "conductorPotentials": [0.9213984270827741, 1.0], + "expansionCenter": [0.0023197106176927533, -0.004847736541849076], + "innerRegionAveragePotential": 0.7447102973543827, + }, + ], + }, + } + mat_assoc = [{"elementIds": [16, 30], "materialId": 1}] + + report, warn = compute_material_rows(material, mat_assoc, 16) + self.assertIsNone(warn) + self.assertIsNotNone(report) + + rows = {j: (c_val, l_val) for j, c_val, l_val in report["rows"]} + self.assertIn(16, rows) + self.assertIn(30, rows) + + self.assertAlmostEqual(rows[16][0], 1.9116497250688996e-11, delta=1e-20) + self.assertAlmostEqual(rows[16][1], 5.820365736777192e-07, delta=1e-15) + self.assertAlmostEqual(rows[30][0], 5.165814539740492e-11, delta=1e-20) + self.assertAlmostEqual(rows[30][1], 2.1538714707844526e-07, delta=1e-15) + + def test_parse_args_requires_positionals(self): + with contextlib.redirect_stderr(io.StringIO()): + with self.assertRaises(SystemExit): + parse_args([]) + + def test_parse_args_happy_path(self): + args = parse_args(["case.tulip.out.json", "16"]) + self.assertEqual(args.json_path, "case.tulip.out.json") + self.assertEqual(args.i_ref, 16) + + +if __name__ == "__main__": + sys.exit(run()) diff --git a/scripts/tulip_paraview_macro.py b/scripts/tulip_paraview_macro.py new file mode 100644 index 0000000..8e37cc2 --- /dev/null +++ b/scripts/tulip_paraview_macro.py @@ -0,0 +1,223 @@ +""" +Tulip ParaView Macro +==================== +Loads all solved conductor solutions for a Tulip case and creates combined +potential views for each multipolar expansion solution in the .tulip.out.json. + +Usage: + 1. Set CASE_DIR to your case directory. + 2. Run as a ParaView macro (Macros > Add New Macro or run from Python Shell). + +The macro creates: + - Layout "Individual Conductors": all PVD readers for electrostatic and + magnetostatic solved conductors (toggle visibility per conductor). + - For each material with a multipolarExpansion, a layout + "Material___Modes" containing two views: + * Electric modes: combined Phi field weighted by conductorPotentials + for every electric solution; stored as arrays ElectricMode_000, ... + * Magnetic modes: same for magnetic solutions (MagneticMode_000, ...). + Switch between modes in the "Color By" selector in the Properties panel. +""" + +import os +import re +import json + +import paraview.simple as pvs + +# ============================================================================== +# CONFIGURATION +# ============================================================================== +# Assumes that results are in this case folder. +CASE_DIR = "PATH_TO_CASE" + +# ============================================================================== +# Helpers +# ============================================================================== + +def _find_pvd(pv_dir: str, conductor_idx: int, kind: str) -> str | None: + name = f"Conductor_{conductor_idx}_{kind}" + pvd = os.path.join(pv_dir, name, f"{name}.pvd") + return pvd if os.path.exists(pvd) else None + + +def _conductor_index(folder_name: str) -> int | None: + m = re.match(r"Conductor_(\d+)_(electrostatic|magnetostatic)$", folder_name) + return int(m.group(1)) if m else None + + +def _make_combined_filter(readers: list, weights_all: list[list[float]], prefix: str): + """ + Return a ProgrammableFilter that reads Phi from each input, computes + every weighted sum listed in weights_all, and stores results as + _000, _001, ... point arrays. + + readers – list of ParaView source objects (same mesh topology) + weights_all – list of weight vectors, one per solution + prefix – array name prefix, e.g. 'ElectricMode' + """ + weights_repr = repr(weights_all) + + script = f""" +import numpy as np +from vtk.util import numpy_support + +weights_all = {weights_repr} + +output.ShallowCopy(inputs[0].VTKObject) + +for sol_idx, weights in enumerate(weights_all): + result = np.zeros(inputs[0].GetNumberOfPoints()) + for w, inp in zip(weights, inputs): + phi_vtk = inp.PointData.GetArray('Phi') + if phi_vtk is not None: + result += w * numpy_support.vtk_to_numpy(phi_vtk) + arr = numpy_support.numpy_to_vtk(result, deep=True) + arr.SetName(f'{prefix}_{{sol_idx:03d}}') + output.GetPointData().AddArray(arr) +""" + pf = pvs.ProgrammableFilter(Input=readers) + pf.Script = script + pf.RequestInformationScript = "" + pf.RequestUpdateExtentScript = "" + return pf + + +# ============================================================================== +# Main +# ============================================================================== + +pvs._DisableFirstRenderCameraReset() + +case_name = os.path.basename(CASE_DIR) +json_path = os.path.join(CASE_DIR, f"{case_name}.tulip.out.json") +pv_dir = os.path.join(CASE_DIR, "ParaView") + +with open(json_path) as fh: + out_data = json.load(fh) + +materials = out_data.get("materials", []) +mat_assoc = out_data.get("materialAssociations", []) + +# Collect sorted conductor indices present in the ParaView folder +all_indices = sorted( + {_conductor_index(d) for d in os.listdir(pv_dir) if _conductor_index(d) is not None} +) + +# -------------------------------------------------------------------------- +# Layout 1 – Individual conductor solutions +# -------------------------------------------------------------------------- +layout_indiv = pvs.CreateLayout("Individual Conductors") +view_indiv = pvs.CreateRenderView() +pvs.AssignViewToLayout(view=view_indiv, layout=layout_indiv) + +elec_readers: dict[int, pvs.PVDReader] = {} +mag_readers: dict[int, pvs.PVDReader] = {} + +for idx in all_indices: + pvd_e = _find_pvd(pv_dir, idx, "electrostatic") + if pvd_e: + r = pvs.PVDReader(FileName=pvd_e) + pvs.RenameSource(f"Conductor_{idx}_electrostatic", r) + elec_readers[idx] = r + d = pvs.Show(r, view_indiv) + d.Visibility = 0 # hidden by default; user toggles in pipeline browser + + pvd_m = _find_pvd(pv_dir, idx, "magnetostatic") + if pvd_m: + r = pvs.PVDReader(FileName=pvd_m) + pvs.RenameSource(f"Conductor_{idx}_magnetostatic", r) + mag_readers[idx] = r + d = pvs.Show(r, view_indiv) + d.Visibility = 0 + +# Show the first electrostatic conductor as a quick default +if all_indices and all_indices[0] in elec_readers: + first_r = elec_readers[all_indices[0]] + disp = pvs.Show(first_r, view_indiv) + disp.Visibility = 1 + pvs.ColorBy(disp, ("POINTS", "Phi")) + pvs.UpdateScalarBars() + +pvs.ResetCamera(view_indiv) + +# -------------------------------------------------------------------------- +# For each material – build combined multipolar views +# -------------------------------------------------------------------------- +for material in materials: + mat_id = material["id"] + mat_type = material.get("type", f"material_{mat_id}") + mp = material.get("multipolarExpansion", {}) + + assoc = next((a for a in mat_assoc if a["materialId"] == mat_id), None) + if assoc is None: + print(f" [warn] No materialAssociation found for material id={mat_id}; skipping.") + continue + + element_ids: list[int] = assoc["elementIds"] + + e_readers = [elec_readers[i] for i in element_ids if i in elec_readers] + m_readers = [mag_readers[i] for i in element_ids if i in mag_readers] + + electric_solutions = mp.get("electric", []) + magnetic_solutions = mp.get("magnetic", []) + + if not electric_solutions and not magnetic_solutions: + continue + + layout_name = f"Material_{mat_id}_{mat_type}_Modes" + layout_mp = pvs.CreateLayout(layout_name) + + # ---- Electric combined view ------------------------------------------ + if electric_solutions and e_readers: + e_weights = [sol["conductorPotentials"] for sol in electric_solutions] + + e_filter = _make_combined_filter(e_readers, e_weights, "ElectricMode") + pvs.RenameSource(f"Material_{mat_id}_{mat_type}_electric_combined", e_filter) + + view_e = pvs.CreateRenderView() + pvs.AssignViewToLayout(view=view_e, layout=layout_mp) + view_e.ViewSize = [800, 600] + + disp_e = pvs.Show(e_filter, view_e) + # Default: show first mode + pvs.ColorBy(disp_e, ("POINTS", "ElectricMode_000")) + pvs.UpdateScalarBars() + pvs.ResetCamera(view_e) + + print( + f" [material {mat_id}] Created electric combined view with " + f"{len(electric_solutions)} modes for {len(e_readers)} conductors." + ) + + # ---- Magnetic combined view ------------------------------------------ + if magnetic_solutions and m_readers: + m_weights = [sol["conductorPotentials"] for sol in magnetic_solutions] + + m_filter = _make_combined_filter(m_readers, m_weights, "MagneticMode") + pvs.RenameSource(f"Material_{mat_id}_{mat_type}_magnetic_combined", m_filter) + + # Split the layout horizontally and add a second view + if electric_solutions and e_readers: + pvs.SplitViewHorizontal(view=view_e, layout=layout_mp) + view_m = pvs.CreateRenderView() + pvs.AssignViewToLayout(view=view_m, layout=layout_mp, hint=2) + else: + view_m = pvs.CreateRenderView() + pvs.AssignViewToLayout(view=view_m, layout=layout_mp) + + view_m.ViewSize = [800, 600] + disp_m = pvs.Show(m_filter, view_m) + pvs.ColorBy(disp_m, ("POINTS", "MagneticMode_000")) + pvs.UpdateScalarBars() + pvs.ResetCamera(view_m) + + print( + f" [material {mat_id}] Created magnetic combined view with " + f"{len(magnetic_solutions)} modes for {len(m_readers)} conductors." + ) + +pvs.Render() +print("Tulip ParaView macro completed successfully.") +print("Tip: In the 'Modes' layout, use 'Color By' to switch between ElectricMode_000, " + "ElectricMode_001, ... arrays to inspect individual multipolar solutions.") diff --git a/src/Launcher.cpp b/src/Launcher.cpp index 8960892..72c230a 100644 --- a/src/Launcher.cpp +++ b/src/Launcher.cpp @@ -63,12 +63,12 @@ Launcher::Launcher(const std::string& inputFile, const std::string& exportFolder void Launcher::run() { std::cout << "Loading input file: " << inputFile_ << std::endl; - const std::string outputPrefix = extractCaseName(inputFile_) + "."; - const std::string driverExportFolder = ensureTrailingSlash(exportFolder_) + outputPrefix; + const std::string caseNamePrefix = extractCaseName(inputFile_) + "."; + const std::string outputPathPrefix = ensureTrailingSlash(exportFolder_) + caseNamePrefix; if (isAdaptedJson(inputFile_)) { auto driver = Driver::loadFromAdaptedFile(inputFile_); - driver.setExportFolder(driverExportFolder); + driver.setExportFolder(outputPathPrefix); std::cout << "Running Tulip analysis..." << std::endl; driver.run(); } @@ -82,7 +82,7 @@ void Launcher::run() Adapter adapter(inputFile_); AdaptedInputParser parser(inputFile_, adapter.getAdaptedInputJSON()); Driver driver(parser.readModel(), parser.readDriverOptions()); - driver.setExportFolder(driverExportFolder); + driver.setExportFolder(outputPathPrefix); std::cout << "Running Tulip analysis..." << std::endl; driver.run(); } diff --git a/src/adapter/Adapter.cpp b/src/adapter/Adapter.cpp index abf1571..b57108c 100644 --- a/src/adapter/Adapter.cpp +++ b/src/adapter/Adapter.cpp @@ -566,6 +566,38 @@ void validateLayerMaterialIds(const nlohmann::json& inputJson) } } +void validateRequiredInputSections(const nlohmann::json& inputJson) +{ + const bool hasMaterials = + inputJson.contains("materials") && inputJson["materials"].is_array(); + const bool hasLayers = inputJson.contains("layers") && inputJson["layers"].is_array(); + + if (hasMaterials && hasLayers) { + return; + } + + std::string missingSections; + if (!hasMaterials) { + missingSections += "'materials'"; + } + if (!hasLayers) { + if (!missingSections.empty()) { + missingSections += " and "; + } + missingSections += "'layers'"; + } + + std::string message = + "Invalid input JSON: missing required top-level array section(s): " + + missingSections + "."; + if (inputJson.contains("model")) { + message += " Found 'model'; expected top-level 'materials' and 'layers'."; + } else { + message += " Expected top-level 'materials' and 'layers'."; + } + throw std::runtime_error(message); +} + std::vector buildAcceptedStepNamesForLayer( const nlohmann::json& layer, const std::map& materialTypeById) @@ -799,6 +831,7 @@ void Adapter::initialize(const nlohmann::json& inputJson, caseName_ = caseName; inputDir_ = inputDir; + validateRequiredInputSections(inputJson); validateLayerMaterialIds(inputJson); adapterOptions_ = parseAdapterOptions(inputJson, std::filesystem::path(inputDir_), caseName_); diff --git a/src/adapter/ShapesClassification.cpp b/src/adapter/ShapesClassification.cpp index 42d7464..dd11ed6 100644 --- a/src/adapter/ShapesClassification.cpp +++ b/src/adapter/ShapesClassification.cpp @@ -334,6 +334,7 @@ bool ShapesClassification::isOpenProblem() const auto roots = nestedGraph.roots(); if (open.size() == 1) return true; + if (roots.empty()) return true; if (roots.size() > 1) return true; if (!roots.empty()) { const auto& root = roots[0]; @@ -415,7 +416,15 @@ EntityMap ShapesClassification::buildVacuumDomain() { EntityMap ShapesClassification::buildClosedVacuumDomain() { const auto roots = nestedGraph.roots(); + if (roots.empty()) { + throw std::runtime_error( + "Unable to build closed vacuum domain: no root entity found."); + } const auto& root = roots[0]; + if (!conductors.count(root)) { + throw std::runtime_error( + "Unable to build closed vacuum domain: root entity is not a conductor."); + } EntityList dom = conductors.at(root); EntityList toRemove; diff --git a/test/adapter/AdapterTest.cpp b/test/adapter/AdapterTest.cpp index 9b5f00a..0fac84a 100644 --- a/test/adapter/AdapterTest.cpp +++ b/test/adapter/AdapterTest.cpp @@ -222,6 +222,39 @@ TEST_F(AdapterTest, dielectric_unshielded_pair_fails_if_step_layer_is_not_presen std::runtime_error); } +TEST_F(AdapterTest, two_wires_open_fails_if_layers_section_is_missing) +{ + const std::string caseName = "two_wires_open"; + nlohmann::json inputJson = readInputJsonFromCaseName(caseName); + inputJson.erase("layers"); + inputJson["model"] = {{"materials", nlohmann::json::array()}}; + + try { + Adapter adapter(inputJson, caseName, inputFolderFromCaseName(caseName)); + FAIL() << "Expected runtime_error"; + } catch (const std::runtime_error& err) { + const std::string message = err.what(); + EXPECT_NE(message.find("'layers'"), std::string::npos); + EXPECT_NE(message.find("top-level 'materials' and 'layers'"), std::string::npos); + } +} + +TEST_F(AdapterTest, two_wires_open_fails_if_materials_section_is_missing) +{ + const std::string caseName = "two_wires_open"; + nlohmann::json inputJson = readInputJsonFromCaseName(caseName); + inputJson.erase("materials"); + + try { + Adapter adapter(inputJson, caseName, inputFolderFromCaseName(caseName)); + FAIL() << "Expected runtime_error"; + } catch (const std::runtime_error& err) { + const std::string message = err.what(); + EXPECT_NE(message.find("'materials'"), std::string::npos); + EXPECT_NE(message.find("top-level 'materials' and 'layers'"), std::string::npos); + } +} + TEST_F(AdapterTest, dielectric_unshielded_pair) { const std::string caseName = "dielectric_unshielded_pair"; @@ -395,3 +428,20 @@ TEST_F(AdapterTest, overlapping_dielectrics_prioritize_higher_relative_permittiv EXPECT_NEAR(rightHighLeftMass, 4.0, 1e-9); EXPECT_NEAR(rightHighRightMass, 8.0, 1e-9); } + +TEST_F(AdapterTest, shapes_classification_without_roots_is_treated_as_open_problem) +{ + gmsh::clear(); + gmsh::model::add("no_roots_case"); + + const EntityList shapes = {}; + const nlohmann::json inputJson = { + {"materials", nlohmann::json::array()}, + {"layers", nlohmann::json::array()} + }; + + ShapesClassification classification(shapes, inputJson); + + EXPECT_TRUE(classification.isOpenCase); + EXPECT_TRUE(classification.isOpenProblem()); +} diff --git a/testData/agrawal1981/agrawal1981.msh b/testData/agrawal1981/agrawal1981.msh index d49f2e7..9d7b1e8 100644 --- a/testData/agrawal1981/agrawal1981.msh +++ b/testData/agrawal1981/agrawal1981.msh @@ -29613,7 +29613,7 @@ $Nodes 29597 -0.8995522719202876 1.763043536181375 0 29598 -0.8877177445822432 1.769014368094762 0 29599 -0.8854043705070319 1.781746044666594 0 -29600 -0.8972652804878468 1.775827504646807 0 +29600 -0.8972652804878468 1.775827504646806 0 29601 -1.310615332880854 0.9476094369748136 0 29602 -1.297944068191723 0.9458983500999457 0 29603 -1.307559396848753 0.9695504897111458 0