segment_network.cpp

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Sebastian Reiter
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */


#include "segment_network.h"
#include "lib_grid/selector.h"
#include "lib_grid/algorithms/geom_obj_util/edge_util.h"
#include "lib_grid/refinement/regular_refinement.h"
#include "lib_grid/file_io/file_io.h"
#include "lib_grid/iterators/lg_for_each.h"

namespace ug{
namespace d3f{
	
/**
 * The last vertex coordinate is interpreted as gauge at each vertex.
 * Upon load, the network is refined until all edges are shorter than maxEdgeLen.
 * Note that the grid 3d coordinates are converted to the 2d (x,y)-coordinates
 * and the gauge.
 */
void SegmentNetwork::load_from_file
(
	const char* name, ///< file name
	number maxEdgeLen ///< prescribed maximum edge length (the grid is refined until this is reached)
)
{
//	use a temporary position attachment to import coordinates and gauge (stored
//	in the last coordinate of each vertex)
	APosition aCoord(aPosition);
	m_grid.attach_to_vertices(aCoord);
	Grid::VertexAttachmentAccessor<APosition> aaCoord(m_grid, aCoord);
	if(!LoadGridFromFile(m_grid, m_SH, name, aCoord)){
		m_grid.detach_from_vertices(aCoord);
		UG_THROW("Couldn't read segment data from '" << name << "'");
	}

//	repeatedly refine all edges which are longer than maxEdgeLen
	if(maxEdgeLen > 0){
		number maxEdgeLenSq = sq(maxEdgeLen);

		Selector sel(m_grid);

	//	helper for more efficient refinement
		AInt aInt;
		m_grid.attach_to_edges(aInt);

		while(1){
			sel.clear();
			lg_for_each(Edge, e, m_grid){
				if(EdgeLengthSq(e, aaCoord) > maxEdgeLenSq)
					sel.select(e);
			}lg_end_for;

			if(sel.empty<Edge>())
				break;
			else
				Refine(m_grid, sel, aInt);
		}

		m_grid.detach_from_edges(aInt);
	}

//	extract position and gauge from aCoord
	lg_for_each(Vertex, v, m_grid){
		m_aaPos[v].x() = aaCoord[v].x();
		m_aaPos[v].y() = aaCoord[v].y();
		m_aaGauge[v] = aaCoord[v].z();
	}lg_end_for;

	m_grid.detach_from_vertices(aCoord);
}

/**
 * Sets the gauge values from a UserData object to the vertices of the segment grid.
 * Note that the userdata object gets the old gauge value as the z coordinate and the
 * subset id of the vertex.
 */
void SegmentNetwork::set_gauge_values
(
	SmartPtr<UserData<number, dim> > h ///< the data to write to the gauge at the vertices
)
{
	MathVector<dim> x;
	
//	loop the vertices
	lg_for_each(Vertex, v, m_grid){
		number& g = m_aaGauge[v];
		x[0] = m_aaPos[v].x(); x[1] = m_aaPos[v].y(); x[2] = g;
		(*h) (g, x, 0, m_SH.get_subset_index(v));
	}lg_end_for;
}

/**
 * Saves the segment grid to a file. The coordinates of the vertices are converted to
 * (x,y,gauge) 3d coordinates.
 */
void SegmentNetwork::save_to_file_with_
(
	gauge_accessor_t& z_acc, ///< to save in the z-coordinate
	const char* name ///< output file name
)
{
	APosition aCoord(aPosition);
	m_grid.attach_to_vertices(aCoord);
	Grid::VertexAttachmentAccessor<APosition> aaCoord(m_grid, aCoord);
	
//	extract position and gauge from aCoord
	lg_for_each(Vertex, v, m_grid){
		aaCoord[v].x() = m_aaPos[v].x();
		aaCoord[v].y() = m_aaPos[v].y();
		aaCoord[v].z() = z_acc[v];
	}lg_end_for;

	SaveGridToFile(m_grid, name, aCoord);

	m_grid.detach_from_vertices(aCoord);
}

/**
 * Finds the vertices in the specified sphere to vrtsOut
 * @return	returns the number of vertices found.
 */
size_t SegmentNetwork::find_vertices
(
	std::vector<Vertex*>& vrtsOut,
	const vector_t& center,
	number radius
) const
{
//	todo: this should be tuned using a quadtree
	number radSq = sq(radius);
	vrtsOut.clear();
	lg_for_each_const(Vertex, v, m_grid){
		if(VecDistanceSq(m_aaPos[v], center) <= radSq)
			vrtsOut.push_back(v);
	}lg_end_for;
	return vrtsOut.size();
}

}//	end of namespace	
}//	end of namespace