errorEstimator delta;
deltaField water.alpha;

// Loops over mesh faces, computing the difference of water.alpha
// between sides of the face
// but caps the difference to 0.5; <== this is hard-coded for now
// on each face, e = max(e, min(mag(ownValue - neiValue), 0.5)) is the error value picked up
// to be normalized by the following settings
// What's important to understand:
// - delta will prefer to refine (a little greedy)
// - it handles boundary patches similarly

// refine if e is between
lowerRefineLevel 0.01;
upperRefineLevel 0.99;
// unrefine if e is less than
unrefineLevel   0.01;
// unrefine if e is greater than
upperUnrefineLevel 0.99;
// else, the cell is queued to retain its current level of refinement
// but this may not be honored due to buffer layers, or different protection settings

errorEstimator  coded;
name emulateDelta; // <== arbitrary name for recognition in dynamicCode folder

// Provide the code to select cells for refinement/unrefinement operations.
// - This only "hints" to the refinement engine what cells to refine/unrefine, it's up
//   to the engine to honor these hints as there a few rules it must follow to produce
//   a valid mesh
// - Other OpenFOAM codeStream constructs (eg. codeInclude, codeOptions, codeLibs)
//   are naturally allowed so you can include headers and link against external
//   libraries freely
code
#{
    Info<< "---->! custom error estimator !<----" << endl;
    error_ = 0.0; // <== set default behavior to not touch the cells
    /*
        End result must be (error_ is a volScalarField here):
        - error_ == 1  if the cell needs to be refined
        - error_ == 0  if the cell is to be left alone
        - error_ == -1 if the vertices of the cell are to be unrefined
    */

    // get references to few useful things
    const auto& alpha = mesh_.lookupObject<volScalarField>("alpha.water");
    const auto& owner = mesh_.owner();
    const auto& neighbour = mesh_.neighbour();
    const label nInternalFaces = mesh_.nInternalFaces();
    // Loop through interal faces and store difference between owner and neighbour
    // note that we must take of the current error value as a single cell will have
    // multiple faces, hence, will be processed multiple times
    for (label facei = 0; facei < nInternalFaces; facei++)
    {
        label own = owner[facei];
        label nei = neighbour[facei];
        scalar eT = mag(alpha[own] - alpha[nei]);
        error_[own] = max(error_[own], eT);
        error_[nei] = max(error_[nei], eT);
    }

    // Settings for "refining everything where delta > 0.01, unrefine otherwise"
    // You only need to set these if you will call the normalize method
    // otherwise, mapping error_ to -1, 0, and 1 manually is perfectly fine
    lowerRefine_ = 0.01;
    upperRefine_ =  GREAT;
    lowerUnrefine_ = lowerRefine_;
    upperUnrefine_ =  GREAT;
    // This converts error_ to the right values (-1, 0, 1) based on the
    // settings lowerRefine_ ... etc and takes care of refining the probes
    if (scale) normalize(error_);
    // error_.correctBoundaryConditions() will be called after this code in case you
    // need to manipulate boundary values
    Info<< "---->! end    error estimator !<----" << endl;
#};

// Include necessary headers, will make included entities available everywhere else
codeInclude #{
    #include "HashTable.H"
#};

// Define custom member data, this is a good place to put variables that
// would live for the whole simulation
// Note: currently, the data defined here must be null-constructed
codeData #{
    HashTable<label> ht_;
#};

// The following code blocks can use most of refinement settings
// mesh_: the mesh object, as const reference to fvMesh
// dict_: a copy of dynamicMeshDict
// error_: the indicator field as a volScalarField
// maxLevel_: the max refinement level, dynamically changeable
// minDx_: the minimal cell size to reach 
// maxRefinementLevel_: the initial max refinement level, a const
// protectedPatches_: as a wordList of patches to protect from refinement
// nPatchesBuffers_: number of cells to protect away from the selected patches
// lowerRefine_, upperRefine_, lowerUnrefine_ and upperUnrefine_

// unless it's annotated as a "const", the user is free to change the values
// of these variables dynamically. The changes will take effect immediately.

// Read-in data from the dynamicMeshDict dictionary if necessary
// the dictionary variable is named "dict"
// this is different from dict_ (which is a copy). "dict", on the other hand
// can pick-up runtime changes to dynamicMeshDict
codeRead #{
    ht_ = dict.lookupOrDefault<HashTable<label>>("theHT", ht_);
#};

// Actual code to set error_ values
code #{
    // use of normalize method is recommended as it offloads the mapping to -1,0,1 values
    // to dictionary settings but this is completely optionally
#};

// You can then put your new "needed" dictionary entries in the dynamicMeshDict

// Eg. Satisfy the ht read
ht ((item1 1) (item2 20));