sohle.h

Go to the documentation of this file.

//
// Created by julian on 9/20/24.
//
 
 
 
#ifndef D3F_UG4_SUBMODULE_SOHLE_H
#define D3F_UG4_SUBMODULE_SOHLE_H
#include <sstream>
 
#include "lib_disc/local_finite_element/lagrange/lagrangep1.h"
#include "lib_disc/spatial_disc/user_data/std_glob_pos_data.h"
namespace ug {namespace d3f{
template <int dim>
class Sohle : public StdGlobPosData<Sohle<dim>, number, dim> {
public:
    virtual void operator() (number& value, const MathVector<dim>& globIP, number time, int si) const {
        evaluate(value, globIP, time, si);
    }
 
    virtual void operator()(number vValue[], const MathVector<dim> vGlobIP[],
                            number time, int si, const size_t nip) const
    {
        for (size_t ip = 0; ip < nip; ++ip) {
            vValue[ip] = vGlobIP[ip][dim - 1]; // Grab the highest dimension coordinate
        }
    }
 
    inline void evaluate(number& value, const MathVector<dim>& globIP, number time, int si) const {
        value = globIP[dim - 1];
    }
 
    template <int refDim>
    inline void evaluate(number vValue[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
        for (size_t ip = 0; ip < nip; ++ip) {
            vValue[ip] = vGlobIP[ip][dim - 1];
        }
    }
};
 
template<typename TDomain>
class ConstantSohlGradient : public StdDataLinker<ConstantSohlGradient<TDomain>, number, TDomain::dim> {
public:
    static const int dim = TDomain::dim;
    typedef typename TDomain::grid_type TGrid;
    typedef typename TGrid::template EdgeAttachmentAccessor<ANumber> accessor_type;
    typedef typename TGrid::template VertexAttachmentAccessor<ANumber> v_accessor_type;
 
    ConstantSohlGradient(const char* subset, number slope, SmartPtr<TDomain> dom) {
        auto grid = dom->grid();
        auto sh = dom->subset_handler();
        auto vAcc = dom->position_accessor();
 
        // 1. Attach the elevation data to the grid vertices
        grid->attach_to_vertices(m_aSohle);
        m_aaSohle.access(*grid, m_aSohle);
 
        // 2. Parse the starting subsets (comma-separated string)
        std::stringstream ss(subset);
        std::string s;
        std::vector<int> subsets;
        while (getline(ss, s, ',')) {
            subsets.push_back(sh->get_subset_index(s.c_str()));
        }
 
        // 3. Initialize BFS structures
        BoolMarker visited(*grid);
        visited.clear();
 
        // Queue stores pairs of: {VertexPointer, AccumulatedDistance}
        std::queue<std::pair<Vertex*, number>> queue;
 
        // Add all vertices in the starting subsets to the queue
        for (int subset_id : subsets) {
            auto vIt = sh->template begin<Vertex>(subset_id, 0);
            auto vEnd = sh->template end<Vertex>(subset_id, 0);
            for (; vIt != vEnd; ++vIt) {
                queue.push({*vIt, 0.0});
                visited.mark(*vIt);
            }
        }
 
        // 4. Perform BFS Traversal
        std::vector<Vertex*> neighborhood;
        while (!queue.empty()) {
            Vertex* v = queue.front().first;
            number curr_dist = queue.front().second;
            queue.pop();
 
            // Set the bed elevation at this vertex
            m_aaSohle[v] = slope * curr_dist;
 
            // Find all connected neighbors
            CollectNeighbors(neighborhood, *grid, v);
 
            for (Vertex* neighbor : neighborhood) {
                if (!visited.is_marked(neighbor)) {
                    visited.mark(neighbor);
 
                    // Get horizontal coordinates (ignore current Z to calculate pure 2D path length)
                    MathVector<dim> v0 = vAcc[v];
                    MathVector<dim> v1 = vAcc[neighbor];
                    v0[dim - 1] = 0.0;
                    v1[dim - 1] = 0.0;
 
                    // Calculate 2D distance and push to queue
                    number step_dist = VecDistance(v0, v1);
                    queue.push({neighbor, curr_dist + step_dist});
                }
            }
        }
    }
 
    inline void evaluate(number& value, const MathVector<dim>& globIP, number time, int si) const {
        UG_THROW("ConstantSohlGradient: Not implemented for single global positions.\n");
    }
 
    template <int refDim>
    inline void evaluate(number vValue[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
        const size_t numSH = 2; // Linear 1D elements have 2 shape functions
        const ReferenceObjectID roid = elem->reference_object_id();
 
        const LocalShapeFunctionSet<refDim>& rTrialSpace =
             LocalFiniteElementProvider::get<refDim>(roid, LFEID(LFEID::LAGRANGE, refDim, 1));
 
        std::vector<number> vShape;
 
 
        Edge* edgeElem = static_cast<Edge*>(elem);
 
        // Loop over integration points
        for (size_t ip = 0; ip < nip; ++ip) {
            rTrialSpace.shapes(vShape, vLocIP[ip]);
            vValue[ip] = 0.0;
 
            // Loop over shape functions to interpolate the vertex data
            for (size_t sh = 0; sh < numSH; ++sh) {
                vValue[ip] += m_aaSohle[edgeElem->vertex(sh)] * vShape[sh];
            }
        }
    }
 
    template <int refDim>
    void eval_and_deriv(number vValue[],
                        const MathVector<dim> vGlobIP[],
                        number time, int si,
                        GridObject* elem,
                        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 = NULL) const
    {
        evaluate<refDim>(vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        // If bDeriv is true, the vvvDeriv matrices remain zeroed out by the caller, which is correct.
    }
 
    bool requires_grid_fct() const { return true; }
 
private:
    ANumber m_aSohle;
    v_accessor_type m_aaSohle;
};
 
 
} // namespace ug
}
#endif //UG4_SOHLE_H