Pe laser scan

.run/pytest test_field.run.xml

@@ -1,13 +1,16 @@
 <component name="ProjectRunConfigurationManager">
   <configuration default="false" name="pytest test_field" type="tests" factoryName="py.test">
     <module name="bgem" />
+    <option name="ENV_FILES" value="" />
     <option name="INTERPRETER_OPTIONS" value="" />
     <option name="PARENT_ENVS" value="true" />
     <option name="SDK_HOME" value="/usr/bin/python3.8" />
     <option name="WORKING_DIRECTORY" value="" />
-    <option name="IS_MODULE_SDK" value="false" />
+    <option name="IS_MODULE_SDK" value="true" />
     <option name="ADD_CONTENT_ROOTS" value="true" />
     <option name="ADD_SOURCE_ROOTS" value="true" />
+    <option name="DEBUG_JUST_MY_CODE" value="true" />
+    <option name="RUN_TOOL" value="" />
     <option name="_new_keywords" value="&quot;&quot;" />
     <option name="_new_parameters" value="&quot;-s&quot;" />
     <option name="_new_additionalArguments" value="&quot;&quot;" />

.run/pytest test_gmsh.run.xml

@@ -1,13 +1,16 @@
 <component name="ProjectRunConfigurationManager">
   <configuration default="false" name="pytest test_gmsh" type="tests" factoryName="py.test">
     <module name="bgem" />
+    <option name="ENV_FILES" value="" />
     <option name="INTERPRETER_OPTIONS" value="" />
     <option name="PARENT_ENVS" value="true" />
     <option name="SDK_HOME" value="/usr/bin/python3.8" />
     <option name="WORKING_DIRECTORY" value="" />
-    <option name="IS_MODULE_SDK" value="false" />
+    <option name="IS_MODULE_SDK" value="true" />
     <option name="ADD_CONTENT_ROOTS" value="true" />
     <option name="ADD_SOURCE_ROOTS" value="true" />
+    <option name="DEBUG_JUST_MY_CODE" value="true" />
+    <option name="RUN_TOOL" value="" />
     <option name="_new_keywords" value="&quot;&quot;" />
     <option name="_new_parameters" value="&quot;&quot;" />
     <option name="_new_additionalArguments" value="&quot;&quot;" />

AGENTS.md

@@ -0,0 +1,96 @@
+# CODEX context for BGEM library
+
+## Project summary
+
+Bgem library is a intended as an extension or wrapper of GMSH python API.
+It helps to keep tracks of the geometry objects during geometrical operations;
+it introduces `ObjectSet` class keeping the tracks of dimtags and provides wrapper
+functions to GMSH API calls.
+
+Further it provides better interface to GMSH fields, reader/writer for GMSH meshes
+and a HealMesh class that can help fixing degenerated elements in a given mesh.
+
+Submodule called `Stochastic` provides generator for DFN (Discrete Fracture Network)
+based on distributions used in geology.
+
+## Project summary (short)
+
+Robust open source tool for creation of parametric geometries and computational meshes via. Python code. 
+Primary focus are hydrogeological applications with geometries including both random fractures and deterministic natural or anthropogenic features.
+Bgem library is a intended as an extension and/or wrapper of GMSH python API.
+
+## CODEX Ignore Folders
+- Not specified
+
+## CODEX Readonly Folders
+- Not specified
+
+## Project Structure
+- `src/bgem` - source files
+- `tests` - pytest based unit tests of individual source blocks
+- `tutorials` - integrated tests, currently insufficient
+
+## `src/bgem/` Project Source Structure
+
+- `bspline/` - 
+- `geometry/` - 
+- `gmsh/` - main functions providing interface for GMSH API
+- `polygons/` - 
+- `stochastic/` - code for DFN generation
+
+
+## CODEX Guidelines
+- treat keyword 'AGENT:' in comments as a source context dependent message for your further development
+- any comment containing `AGENT:` is an active developer instruction
+- NEVER remove, rewrite, or move an `AGENT:` comment unless you implement that instruction in the same change
+- if you make only a local fix around an `AGENT:` comment, leave the comment untouched
+- if an `AGENT:` instruction looks outdated or wrong, ask before removing it
+- never stage or commit changes yourself
+- Always review your changes before finishing for human review.
+- Do not touch anything outside the repository directory, unless directly asked!
+
+## Status tracking
+- `STATUS.md` is the handoff log for interrupted or multi-turn work. Update it when a session ends with unfinished relevant work, when the user asks for a status review, or when a fresh checkpoint would help the next session continue without re-discovery.
+- Keep newest entry first. Do not rewrite older entries except to fix clearly factual mistakes.
+- Start each entry with one line in this form:
+  `` `YYYY-MM-DD`: `<commit>` @ `<branch>` by `<author>` ``
+- Under each entry keep exactly these sections in order:
+  `## Goal`
+  short statement of the intended task or checkpoint scope
+  `## Changes summary`
+  flat bullets describing committed changes first, then important staged/unstaged/untracked changes if they are relevant to continuing the task
+  `## Verified`
+  flat bullets with commands actually run and the important observed result
+  `## Open items`
+  flat bullets for remaining risks, missing verification, known breakage, or next recommended step
+- Record only repo-relevant facts that help continuation. Skip conversational history, speculation, and incidental noise.
+- When the worktree is dirty, explicitly distinguish committed, staged, unstaged, and untracked changes, but mention only files relevant to the tracked task.
+- Prefer clickable file links for important files mentioned in `STATUS.md`.
+- If verification was partial, say so plainly. Do not imply a full test pass when only compile checks, a single testcase, or CLI smoke checks were run.
+- If a run failed, record the failing command and the actionable failure mode instead of hiding it.
+- Before finishing a task that changed the practical project state, review whether `STATUS.md` still matches the actual branch/worktree state and update it if needed.
+
+## Coding rules
+- Best code, is no code!
+- prefer functional style with poor functions;
+  ideally do not change objects after construction, all methods do calculations
+  only reading the data in the class
+- prefer high level code: numpy, pandas, xarray instead loops and native python structures (lists, dicts)
+- use logging
+- Use logging for debug outputs.
+- use pathlib
+- use attrs for dataclasses
+- use attrs staticmethod/classmethod technique to construct from other data then is stored in the dataclass
+- Be defensive, with strong checks, but only for the user input data.
+  That means error inputs must raise early. Therefore only check for existing keys in input dicts
+  if these will be required down in a long calculation. Otherwise just let KeyError do the job.
+- do not use "guess" default values, only obvious defaults
+- Do not use other config keys in the case of a KeyError, just throw early.
+- Do just basic asserts for consistency for function inputs.
+- Can add more asserts if needed during debugging.
+- NEVER resolve test errors by try blocks
+- NEVER add runtime fallbacks or import shims to compensate for a broken or incomplete environment.
+  If a declared dependency or tool is missing, report the environment problem plainly and fix the environment or tests around it, but do not implement code workarounds.
+- NEVER write "self explanatory" into comments
+- in comments indicate by ?? if you are not certain about intent of particular variable, function, parameter ...
+

requirements.txt

@@ -3,3 +3,4 @@ numpy
 scipy
 bih
 plotly
+attrs

setup.py

@@ -56,7 +56,7 @@
     include_package_data=True,
     zip_safe=False,
     #install_requires=['numpy', 'scipy', 'bih', 'gmsh-sdk<=4.5.1'],
-    install_requires=['numpy>=1.13.4', 'pandas', 'scipy', 'bih', 'gmsh>=4.10.4'],
+    install_requires=['numpy>=1.13.4', 'pandas', 'scipy', 'bih', 'gmsh>=4.10.4', 'attrs'],
     # incompatible changes in SDK in release 4.6.0 to be changed in the new release of bgem
     python_requires='>=3',
     # extras_require={

src/bgem/gmsh/field.py

@@ -347,6 +347,21 @@ def construct(self, model):
         # creata MathEval
 
 
+# def constant(value):
+#     """
+#     Define constant field.
+#     Unfortunately much slower then the Box field alternative.
+#     Seems that MathExpr evaluation time is proportional to the string length, i.e. it is parsed every time.
+#     :param value:
+#     :return:
+#     """
+#     return FieldExpr(str(value), [])
+
+def constant(value):
+    """
+    Make a field with constant value = value.
+    """
+    return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value)
 
 """
 Predefined coordinate fields.
@@ -414,11 +429,6 @@ def box(pt_min: Point, pt_max: Point, v_in:float, v_out:float=1e300) -> Field:
 
 
 
-def constant(value):
-    """
-    Make a field with constant value = value.
-    """
-    return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value)
 
 
 
@@ -549,30 +559,41 @@ def minimum(*fields: Field) -> Field:
 # Structured
 
 
-def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
+
+
+def linear(x_field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
     """
-    Geometric interpolation between points a = (x0, y0) and b = (x1, y1).
+    Linear transformation of the scalar field given by two points.
+    Prescribed mapping a[0] to a[1], and b[0] to b[1].
     """
-    x = field
     x0, y0 = a
     x1, y1 = b
-    l0, l1 = math.log(y0), math.log(y1)
-    a = (l1 - l0) / (x1 - x0)
-    x_avg = (l1 * x0 - l0 * x1) / (l1 - l0)
-    return current_module.exp(a * (x - x_avg))
+    a = (y1 - y0) / (x1 - x0)
+    #x_avg = (y1 * x0 - y0 * x1) / (y1 - y0)
+    x_avg = (x0 + x1) / 2.0
+    y_avg = (y0 + y1) / 2.0
+
+    if a > 1e10:
+        # Following formula could provide more stable results
+        # but needs one more scalar addition
+        # that makes the test 3D meshing about 2% slower.
+        return a * (x_field - x_avg) + y_avg
+    else:
+        # Just two field operations.
+        add = - (a * x0 - y0)
+        return a * x_field + add
 
 
-def linear(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
+def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
     """
-    Linear transformation of the scalar field given by two points.
-    Prescribed mapping a[0] to a[1], and b[0] to b[1].
+    Geometric interpolation between points a = (x0, y0) and b = (x1, y1).
+
     """
     x = field
     x0, y0 = a
     x1, y1 = b
-    a = (y1 - y0) / (x1 - x0)
-    x_avg = (y1 * x0 - y0 * x1) / (y1 - y0)
-    return a * (x - x_avg)
+    l0, l1 = math.log(y0), math.log(y1)
+    return current_module.exp(linear(field, (x0, l0), (x1, l1)))
 
 def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: float = 1) -> Field:
     """
@@ -583,9 +604,7 @@ def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: flo
     x0, y0 = a
     x1, y1 = b
     l0, l1 = y0 ** (1/q), y1 ** (1/q)
-    a = (l1 - l0) / (x1 - x0)
-    x_avg = (l1 * x0 - l0 * x1) / (l1 - l0)
-    return (a * (x - x_avg)) ** q
+    return linear(field, (x0, l0), (x1, l1)) ** q
 
 def threshold(field:Field, lower_bound:Tuple[float, float],
               upper_bound:Tuple[float, float]=None, sigmoid:bool=False) -> Field: