This is related to #89, but discusses a different style of refinement, subdividing lines into N segments. This is probably especially useful for Inhoms (e.g. BuildingPits).

The implementation requires some thought, hence this long and detailed post... This is just one idea, but maybe serves as a good starting point for figuring out the best way forward.

Refining elements into segments

The general idea is to give certain elements a Element._refine() method. This method is called in initialize() if refine=True for that particular element. The _refine() method takes n as an argument, which splits line segments into $N$ segments, and returns a list of new smaller (shorter?) elements.

Some requirements:

The refinement is performed internally, the user only needs to indicate it needs to be done with a refine=True/False flag in the element.
Elements without a refine option should continue to work normally.
The user-input to the model should be stored.
Repeated calls to initialize() or solve() should always give the same answer, obviously.

Limitations:

This implementation splits each line segment into $N$ parts, regardless of its length.

Some example code:

ml = tml.ModelMaq(...)
bpit = tml.LeakyBuildingPitMaq(..., refine=True)
well = tml.Well(...)
ml.solve()

ml.elements  # contains Well
ml.aq.inhoms  # contains original LeakyBuildingPitMaq

ml.aq.inhomlist  # contains refined LeakyBuildingPitMaq
ml.elementlist  # all computation elements, so including refined intlinesinks created by inhom

Distinguishing between user-added elements and internally added elements

For this to work, there needs to be a distinction between elements specified by the user and elements that are created internally when refine is called.

Elements defined by user:
- Model.elements includes all elements added to a Model by the user
- Model.aq.inhoms includes all inhoms added to a a Model by the user
Elements used for computation:
- Inhom.elementlist is the list of all elements belonging to a certain inhom for the computation
- Model.elementlist is the list of all elements for the computation
- Model.aq.inhomlist is the list of inhoms used for computation

All elements and inhoms will need an addtomodel option, to prevent internally created elements from getting added to the user element list. Some elements already have this, but now more will need this option.

The different types of elements

Point elements, e.g. wells, are not refined.
Single linesinks must produce new elements with either scaled or copied attributes (e.g. discharge, head, etc.)
PolyLineSinks can be refined by creating a new copy of itself with updated x, y coordinates. This requires moving the code for building the linesinks into a separate function. This function then has to be called post-refinement. These sub-elements are never added to the model, and only exist under the PolyLineSink class.
Inhoms are refined by creating a new copy of itself with updated x, y coordinates. The create_elements method of Inhoms produces a list of elements that are added to the computation element list.

Workflow

For elements

Model.add_element() adds element to Model.elements
- For simple elements this is a straight-forward add
- For compound elements, the compound element is added, not its sub-elements.
- This list tracks all user-added elements.
Model.initialize()
- Model.elementlist is populated in this step, for collecting all computation elements.
- Optionally calls refine. Refine creates new elements (which are not added to the model in constructor, addtomodel=False) that are added to Model.elementlist for the computation.
- Adds original or refined elements to model.aq.elementlist.
Model.elementlist and Model.aq.elementlist start empty at the start of every solve.

For inhoms

Model.aq.add_inhom() adds inhom to the user-inhom list Model.aq.inhoms.
- This provides a list of inhoms added by the user.
Model.aq.initialize() eventually calls inhom.create_elements()
- create_elements creates elements that are added to the computation list Model.elementlist. These elements are added manually, and not by the constructors of the elements themselves (addtomodel=False).
- Optionally calls refine. Refine creates a new refined inhom (this refined inhom is not added to model.aq.inhoms, but added to computation inhom list: model.aq.inhomlist).
- Calls create_elements on refined or original Inhom.

Refinement based on proximity vertices

Link to issue #89.

Additional automatic refinement is desirable when elements lie close to one another, e.g. a well near an impermeable wall. The modifications presented here should also make it possible to implement this style of refinement.

One challenge is that the current segment-splitting refinement will modify the vertices of elements in the computation. These refined elements only become available in the initialize function, so not all vertices may be known yet in the initialize phase. That means a second loop through the element list is required in Model.initialize() to further refine them based on vertices.

# in Model.initialize()

# 1. optionally refine elements and build computation element list

# 2. optionally refine and initialize inhoms, adding inhom elements to computation list

# collect vertices from original elements.
refined_element_list = [] # new list for further refinement
xverts = np.concatenate([e.x for e in self.elements])
yverts = np.concatenate([e.y for e in self.elements])
for e in self.elementlist:
   if hasattr(e, "_refine_vertices"):
      refined_elems = e._refine_vertices(xverts, yverts)
      refined_element_list += refined_elems
   else:
      refined_element_list.append(e)
# overwrite original elementlist
self.elementlist = refined_element_list
for e in self.elementlist:
   e.initialize()

Questions

The user-specified elements can be different from the actual elements used in the computation. Any methods called on the original element will either not work or be incorrect.
- Perhaps we can catch this with the refine flag? If refine=True, the element was not used in the computation, and therefore should throw an error/warning if any methods are called on it?
- some code to illustrate:
```
ml = tml.ModelMaq(...)
hls = tml.HeadLineSink(..., refine=True)
ml.solve()
```
hls.discharge() # <-- this is not the element that was used in the computation! ml.elementlist[0].discharge() # <-- this would be the refined element that was used
Should there be a model-level refine option, e.g. ml.solve(refine=True)? If so, would this approach make sense:
- refine=None, use Element level options (default)
- refine=False, override Element level options with False
- refine=True, override Element level optinos with True
- refine=<int>, add this as an option to set refinement level globally?
Some elements had an if self.addtomodel: self.aq.add_element(self) in the initialize function. Is there ever a reason not to add an element to the aquifer computation list in initialize?
An alternative to creating a new refined inhomogeneity or compound element instance (e.g. HeadLineSinkString) is to refine the element in-place, prior to creating the sub-elements. The advantage is that the original reference to the object can be used for computation, as it will be used in the calculation. The downside is that you lose a copy of the original user-specified element.

mbakker7 / timml

Refining elements into N segments #103