richards_fs_height_impl.h

Go to the documentation of this file.

//
// Created by julian on 2/27/26.
//
 
#ifndef UG4_RICHARDS_FS_HEIGHT_IMPL_H
#define UG4_RICHARDS_FS_HEIGHT_IMPL_H
#include "lib_disc/local_finite_element/local_finite_element_provider.h"
#ifdef UG_PARALLEL
     #include "lib_grid/parallelization/util/compol_attachment_reduce.h"
     #include "lib_grid/parallelization/parallel_grid_layout.h"
#endif
namespace ug {
namespace d3f {
 
// -------------------------------------------------------------------------
// Constructor Implementation
// -------------------------------------------------------------------------
template<typename TGridFunction>
RichardsFreeSurfaceHeight<TGridFunction>::RichardsFreeSurfaceHeight
(
    SmartPtr<gf_type> spGf,
    const char *subset,
    const char *cmp
)
    :m_spGf(spGf),
    m_cmp(cmp),
    m_glob_userdata(spGf, cmp),
    m_sp_ls_pos_z(spGf, spGf->fct_id_by_name(cmp), 0) // Sets up the level set position measurer for z-direction
{
    m_sp_ls_pos_z.reinit();
    // Set up the subset group to restrict operations to the specified boundary (e.g., "top")
    m_SsGrp.set_subset_handler(spGf->domain()->subset_handler());
    m_SsGrp.add(TokenizeString(subset));
    m_sp_vrt_val = make_sp(new std::vector<number>());
    // Attach a vertex numbering property to the grid to uniquely identify vertices during evaluation
    spGf->domain()->grid()->attach_to_vertices_dv(m_a_vrt_num, -1);
    spGf->domain()->grid()->attach_to_vertices_dv(m_a_vrt_val, 0.0);
 
}
 
// -------------------------------------------------------------------------
// Reinit Implementation
// -------------------------------------------------------------------------
    template<typename TGridFunction>
void RichardsFreeSurfaceHeight<TGridFunction>::reinit() {
 
    m_sp_ls_pos_z.reinit();
 
    Grid& grid = *m_spGf->domain()->grid();
 
    // Attachment access for global index and vertex position
    num_aa_type aa_vrt_num (grid, m_a_vrt_num);
    val_aa_type aa_vrt_val (grid, m_a_vrt_val);
    auto vAccPos = m_spGf->domain()->position_accessor();
 
    // Clear previous cached measurements
    m_measurement_points.clear();
    m_sp_vrt_val->clear();
 
    size_t n = 0; // Global index counter for the vertices in the subset
 
    // Vector to collect local vertices to update cached vectors after communication
    std::vector<Vertex*> local_vertices;
 
    // --- Phase 1: Local Setup & Evaluation ---
    for (size_t i = 0; i < m_SsGrp.size (); i++)
    {
        int si = m_SsGrp[i];
 
        for (auto iter = m_spGf->template begin<Vertex> (si);
                    iter != m_spGf->template end<Vertex> (si); ++iter, ++n) {
 
            aa_vrt_num[*iter] = n;
            local_vertices.push_back(*iter);
 
            // Extract the lower-dimensional (x, y) coordinates
            MathVector<dim - 1> xy_coords;
            for (size_t d = 0; d < dim - 1; d++) {
                xy_coords[d] = vAccPos[*iter][d];
            }
 
            // Create a new measurement point and add it to our tracking vector
            m_measurement_points.push_back(xy_coords);
            m_sp_vrt_val->push_back(0); // pre-allocate
 
            auto& pnt = m_measurement_points.back();
            pnt.valid = m_sp_ls_pos_z.get_height_at(pnt.xy, pnt.z);
 
            if (pnt.valid) {
                aa_vrt_val[*iter] = pnt.z;
            } else {
                // If invalid locally, assign a highly negative sentinel for the MAX reduction
                aa_vrt_val[*iter] = -std::numeric_limits<number>::max();
            }
        }
    }
 
#ifdef UG_PARALLEL
    // --- Phase 2: Interface Direct Communication ---
    // Extract layout types and manager strictly for Vertices
    typedef typename GridLayoutMap::Types<Vertex>::Layout layout_t;
    DistributedGridManager& dgm = *grid.distributed_grid_manager();
    GridLayoutMap& glm = dgm.grid_layout_map();
    pcl::InterfaceCommunicator<layout_t> icom;
 
    // We extract the exact Attachment type defined
    typedef decltype(m_a_vrt_val) a_val_type;
 
    // 1. REDUCE: Send local values from Slaves to Master and take the MAX
    // We use MAX so the valid height beats the -MAX sentinel
    ComPol_AttachmentReduce<layout_t, a_val_type> compolMax(grid, m_a_vrt_val, ComPol_AttachmentReduce<layout_t, a_val_type>::MAX);
    icom.exchange_data(glm, INT_H_SLAVE, INT_H_MASTER, compolMax);
    icom.communicate();
 
    // 2. COPY: Broadcast the winning MAX value from Master back to all Slave nodes
    ComPol_CopyAttachment<layout_t, a_val_type> compolCopy(grid, m_a_vrt_val);
    icom.exchange_data(glm, INT_H_MASTER, INT_H_SLAVE, compolCopy);
    icom.communicate();
 
 
#endif
 
    // --- Phase 3: Finalize Caches ---
    // Now that the attachment has been synced natively, we read it back
    // to populate your structs and m_sp_vrt_val for eval_and_deriv
    for (size_t j = 0; j < local_vertices.size(); j++) {
        Vertex* vrt = local_vertices[j];
        number synced_z = aa_vrt_val[vrt];
        auto& pnt = m_measurement_points[j];
 
        if (synced_z == -std::numeric_limits<number>::max()) {
            // No process found the surface at this node
            pnt.valid = false;
            pnt.z = 0.0;
            aa_vrt_val[vrt] = 0.0; // Clean up the sentinel
        } else {
            // Surface was found!
            pnt.valid = true;
            pnt.z = synced_z;
            (*m_sp_vrt_val)[j] = synced_z;
        }
    }
}
 
 
// -------------------------------------------------------------------------
// Evaluation Implementation
// -------------------------------------------------------------------------
template<typename TGridFunction>
template<int refDim>
void RichardsFreeSurfaceHeight<TGridFunction>::eval_and_deriv
(
    number vValue[], 
   const MathVector<dim> vGlobIP[],
   number time,
   int si,
   GridObject* grid_obj,
   const MathVector<dim> vCornerCoords[],
   const MathVector<refDim> vLocIP[],
   const size_t nip,
   LocalVector* u,
   bool bDeriv,
   int s,
   std::vector<std::vector<number> > vvvDeriv[],
   const MathMatrix<refDim, dim>* vJT
)
{
    typedef typename gf_type::template dim_traits<refDim>::grid_base_object elem_type;
 
    // Guard against evaluating on a subset we didn't cache data for
    if (! m_SsGrp.contains (si))
        UG_THROW ("RichardsFreeSurfaceHeight: Attempt to get the value for a subset that has not been specified.");
 
    // Cast to the correct element type and fetch local shape functions for interpolation
    const ReferenceObjectID roid = grid_obj->reference_object_id ();
    const elem_type * elem = dynamic_cast<elem_type *> (grid_obj);
    const LocalShapeFunctionSet<refDim> & rTrialSpace =
        LocalFiniteElementProvider::get<refDim> (roid, LFEID (LFEID::LAGRANGE, refDim, 1));
 
    num_aa_type aa_vrt_num (* m_spGf->domain()->grid (), m_a_vrt_num);
    val_aa_type aa_vrt_val (* m_spGf->domain()->grid (), m_a_vrt_val);
    std::vector<number> shape;
 
    // Loop over all integration points (IPs) within the element
    for (size_t ip = 0; ip < nip; ++ip) {
        // Evaluate the shape functions at the local coordinate of the current IP
        rTrialSpace.shapes (shape, vLocIP [ip]);
 
        // Initialize the interpolated value
        vValue [ip] = 0;
 
        // Interpolate the height using the vertex values and shape functions
        for (size_t sh = 0; sh < shape.size (); ++sh) {
            const Vertex* vrt = elem->vertex (sh);
            size_t vrt_no = aa_vrt_num [vrt]; // Retrieve our cached vertex index
 
            // Check if a valid free surface was found at this vertex during reinit()
            if (m_measurement_points[vrt_no].valid) {
                // Surface found: use the cached free surface height (z)
                //UG_LOG("valid vertex"<<vrt_no<<": " << aa_vrt_val[vrt]<< "\n")
                vValue [ip] += aa_vrt_val[vrt] * shape [sh];
            }
            else {
                // Surface NOT found (e.g., domain is fully saturated and free surface is above it)
                // Fallback: use the physical z-coordinate of the domain boundary vertex
                //UG_LOG("invalid vertex"<<vrt_no<<": " << vCornerCoords[sh][dim -1]<< "\n")
                vValue [ip] += vCornerCoords[sh][dim -1] * shape [sh];
            }
        }
    }
}
}
}
 
 
#endif //UG4_RICHARDS_FS_HEIGHT_IMPL_H